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PREFACE 


The current Final Report contains results of the study which 
was performed in Scientific Research Center "ECOLEN" (Moscow, 
Russia) according to National Aeronautics and Space Administration 
Cooperative Agreement No.NCCW-75. The study was addressed to the 
development and verification of non-expensive approach for 
modeling of supersonic turbulent diffusion flames based on 
flamelet consideration of the chemistry/turbulence interaction (FL 
approach) . Research work included: development of the approach and 
CFD tests of the flamelet model for supersonic jet flames; 
development of the simplified procedure for solution of the 
flamelet equations based on partial equilibrium chemistry 
assumption; study of the flame ignition/extinction predictions 
provided by flamelet model. The performed investigation 
demonstrated that FL approach allowed to describe satisfactory 
main features of supersonic H 2 /air jet flames. Model demonstrated 
also high capabilities for reduction of the computational expenses 
in CFD modeling of the supersonic flames taking into account 
detailed oxidation chemistry. However, some disadvantages and 
restrictions of the existing version of approach were found in 
this study. They were: i) inaccuracy in predictions of the passive 
scalar statistics by our turbulence model for one of the 
considered test cases,- ii) applicability of the available version 
of the flamelet model to flames without large ignition delay 
distance only. 

Based on the results of the performed investigation, we 
formulated and submitted to the National Aeronautics and Space 
Administration our Project Proposal for the next step research 
directed toward further improvement of the FL approach. 
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INTRODUCTION 


Accurate prediction of the turbulent combustible flows 
requires to utilize joint probability density function (pdf) for 
averaging of highly nonlinear chemical source terms in reactive 
scalars conservation equations. Neglecting of the 
turbulence/chemistry interaction effects (so-called "quasilaminar" 
approach) can lead to sufficient inaccuracy in predictions [ 1 , 2 ]. 

Usually, modelers utilize approach where pdf form is assumed 
based on some kind of intuitive consideration (Assumed PDF 
approach) [3-6] . Here, successful choice of the pdf form is based 
fully on the intuitia of modeler and it can not be unique for 
different reacting systems. Much more elaborate way for pdf 
construction is the solution of evolution equation for pdf using 
Monte-Carlo simulation (Evolved PDF approach) [7,8]. Significant 
progress has been achieved during last years due to both computers 
and appropriate numerical algorithms fast improvement [9-11] . 
However, up to now, Evolved PDF modeling requires enormous 
computational expenses due to large multidimensionality of pdf 
evolution equation [12] . 

One of the ways for development of the computationally 
non-expensive procedure for pdf construction is flamelet approach 
[13-16] . The simplification of the problem is achieved here based 
on physical assumption that chemical processes are mostly confined 
in the local vicinity of the near-stoichiometric surfaces (this 
feature is approximately valid for many classes of the turbulent 
flames) . The assumption about small thickness of the reaction 
zones allows to reduce instantaneous mass conservation equations 
for reactive scalars to the system of the ordinary differential 
equations (flamelet equations). Its solution gives relations for 
reactive species mass fractions and temperature depending on 
mixture fraction z and its scalar dissipation N=D(Vz) 2 (D is 
molecular diffusivity) i.e. Ca=Ca(z,N), T=T(z,N). The later 
relations allows to present joint pdf for reactive scalars 
p (C ,...., C ,T) depending on mixture fraction and scalar 
dissipation pdf p(z,N) and to reduce consideration of the reactive 
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scalars statistics to the large- and small-scale statistics of 
passive scalar. The passive scalar pdf is well -investigated and 
it can be modeled based on the first couple of the moments for 
mixture fraction [1,2,17]. 

There is additional advantage of the FL approach. In many 
cases, flamelet equations can be integrated before the start of 
hydrodynamics calculations. So, CFD modeling can be performed 
using tabulated solutions for reactive scalars and temperature 
versus mixture fraction z and scalar dissipation N (flamelet 
library approach) . Thus, a time needed for numerical calculations 
of non-premixed flames could be significantly reduced. By such a 
manner very complex detailed kinetics schemes can be incorporated 
into the CFD codes without computational time increasing. 

The FL approach versions varies between creators [13,18-20] . 
In our studies, we use flamelet model proposed by Dr . V . Kuznetsov 
in [19,21], Previously, Kuznetsov's flamelet model have been 
tested using numerous data obtained in free- and confined subsonic 
jet diffusion flames (H /air, CH /air, C H /air) [1,22]. Flamelet 

2 4 3 8 

model demonstrated quite satisfactory capabilities in predictions 
of temperature, stable species and radicals concentrations. It was 
used also for prediction of nitric oxides (NOx) emissions from 
diffusion flame combustors of gas turbine engines [23] . 

This state of art was the starting point for the current 
one-year investigation. Its goal was to generalize the flamelet 
model approach for new classes of combustible flows and to 
investigate its computational capabilities for CFD modeling of 
diffusion flames. The discussions which we were able to have with 
Dr. L.Povinelli (NASA Lewis Research Center) allowed to adjust the 
problem still further. The supersonic jet configurations were 
chosen for investigation. Test cases [24,25], where supersonic 
H 2 /air jet flames were studied experimentally, were selected for 
model validation. 

Specifically, three tasks were formulated for the current 
investigation : 

o Development of non-expensive approach for modeling of 
supersonic jet flames based on FL model. Verification of its 
computational efficiency and accuracy of predictions in CFD 
tests. 
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« Consideration of the possible simplifications in flamele t 
calculations based on partial equilibrium assumption for the 
detailed oxidation chemistry . 

o Study of the FL predictions features for ignition/extinction 
phenomena in high-enthalpy flows. 

The interim results which were obtained during research work 
implementation were documented in our two Reports to NASA 

[26] , [27] . The main results of the research were reported and 
discussed also during NASA delegation visit to Russia in March 
1996. The current Final Report summarizes results of the performed 
study as a whole. 

The Final Report is organized as follows. 

Part I describes flamelet approach (Sec.I.l), its averaging 
procedure (Sec. I. 2) and developed procedure for the flamelet model 
incorporation into compressible CFD solvers (Sec. 1. 3) . Account of 
the flamelet model equations is given in Appendix A. 

Part II is addressed to the CFD tests of the flamelet model 
capabilities (both computational and physical) . The detailes of 
CFD tests for the conditions of the Beach et al . experiment [24] 
are presented in Sec. II. 1. and for conditions of Burrows -Kurkov 
experiment [25] - in Sec. II. 2. The additional verifications 

performed at the final step of the research does not change 
results and conclusions concerning performed CFD tests reported in 
interim Report [27] . So, Sec. II. 1 and II. 2 compiled mostly results 

[27] . Additional illustrations were introduced in these Sections 
to demonstrate accuracy of performed computations only. Brief 
summary of the CFD tests results, our conclusions concerning model 
capabilities and ways for its improvement are given in Sec. II. 3. 
Comparative estimations of the flamelet model capabilities with 
Evolved and Assumed PDF modeling are presented in Sec. II. 3 also. 
Used computational codes are described in Appendix B. Used 
thermochemistry approximations are presented in Appendix C. 

The Part III contains results of the supplemented studies 
which were obtained at the final step of the research. Sec. III. 1 
describes semi -analytical procedure which was developed for the 
additional reduction of computational expenses and simplification 
of the flamelet calculations based on the partial equilibrium 


6 



assumption for the H 2 oxidation chemistry (two-body reactions are 
equilibrated) . The results of the parametric flamelet model 
calculations of the ignition/extinction phenomena for 

high-enthalpy flames are presented in Sec. III. 2. Conclusions 
concerning results of supplemented studies are given in Sec. III. 3. 

Research team greatly thanks to Dr. Louis Povinelli (NASA 
LeRC) for the formulation of the problem for current investigation 
and for fruitful results discussions during study implementation. 
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Part I. DESCRIPTION OF THE FLAMELET MODEL APPROACH (FL) 


I.l ASSUMPTIONS AND EQUATIONS 

The characteristic feature of the majority of the nonpremixed 
combustion problems is that the maximum temperature and the 
highest reaction rates are observed in the local vicinity of the 
surfaces with the stoichiometric mixture composition. As a rule 
the relative role of the chemical reactions outside 
near- stoichiometric zones is small enough. So one can expect that 
thickness 5 of these near- stoichiometric reactions zones 

ch 

(flamelets) is small enough (see Fig.l). 

In turbulent flows the stoichiometric surface is highly curved 
and randomly fluctuated. To characterize its location the mixture 
fraction z is introduced as the total mass fraction of all kinds 
of atoms initially been contained in fuel and then converted to 
other chemical species arising in the flame. The mixture fraction 
z equals to 0 in the flow of pure oxidizer and it equals to 1 in 
the flow of pure fuel. The mixture fraction z has value 
z=z =l/(l+St) at the stoichiometric surface (dotted line in 

S 

Fig.l), where St is the mass stoichiometric coefficient. Using 
the atoms conservation equations and neglecting the difference in 
molecular dif f usivities of reactive species one has the following 
equation which mixture fraction obeys: 
az 

p-g£- + p(UV)z = VpDVz (I.l) 

-> 

where t is time; p is density; D is molecular diffusivity; U is 
flow velocity. 

The reactive species conservation equations can be seriously 
simplified based on the assumption about small thickness of the 
reaction zones . The convective and unsteady terms can be dropped 
out and mixture fraction z can be used as an independent variable 
instead of space coordinate (full account is somehow length, that 
is why it is given in Appendix A) . As the result, reactive species 
conservation equations are reduced to the following system of 
the ordinary differential equations: 

N s d Ca + Ra = 0 a=l, . . . , J (1.2) 

dz 2 
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where Ca are the reactive species mass fractions; Ra are the 
chemical production terms; J is total number of reactive species; 

s fdz 1 ^ 

parameter N S =D is the value of instantaneous scalar 

dissipation N=D(Vz) 2 at the stoichiometric surface which 
characterize the reactive species fluxes to the reaction zones 
[1,2]; D is molecular diffusivity; n is coordinate normal to the 
surface z=z (Fig.l) . 

S 

The mixture fraction fluctuations have the turbulent integral 
length scaling but the scalar dissipation fluctuations have 
turbulent micro-length scaling [1,2]. So, it is expected that 
mixture fraction and scalar dissipation are non-correlated (inside 
the turbulent mixing layer) and parameter N s is treated in the 
flamelet model equations (1.2) as some random and fluctuating 
number, which does not depend on z. 

The same kind of reasoning can be applied to the energy 
conservation equation. Here the additional suggestions are: 
i) Lewis number Le equal to unity; ii) the role of the unsteady 
pressure fluctuations, viscous dissipation and radiative heat 
losses terms is small enough. As the result the energy 
conservation equation can be reduced to the following form: 

^!h . 0 

J 2 

dz 


(i.3: 


where total enthalpy H is defined as H = h + (U-U)/2; h=^ h C is 

a= i a a 

static enthalpy; species specific enthalpies h^= Jcp^dT + Ah a (To) 

T o 

are used taking into account species heats of formation Ah^ at 
reference temperature T=To . 

The boundary conditions (BC) for the flamelet model equations 
( 1 . 2 ) - ( 1 . 3 ) are posed at z=l (pure fuel) and z=0 (pure oxidizer): 

z=0 H=H a ; C =C A 
' a a 


z=l H=H f ; C =C F 
' a a 


<x=l, 


, J 


(i.4: 


where superscripts f and a denote composition and total 
enthalpies for the flows of fuel and oxidizer respectively. 

The eq. (1.3) is integrated over z from 0 to 1 and the total 
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F A 

enthalpies of hydrogen (H ) at z = l and air (H ) at z=0 is used 
to define the constants in the obtained linear relation. As the 
result the flamelet model equations (1.2), (1.3) are re-written in 
a form: 


N s d Ca + R = o 
dz 2 

h=(H F -H A ) z + H A • 


U 2 /2 


a = l, . . . , J 


(1.5) 


The formulated flamelet model boundary problem 
boundary conditions (1.4) gave the solution for 
temperature T in the following parametric form: 


C =C (z, 
a a ' 


N 


IT 


,BC) 


(1.5) with the 

C and static 

a 


n 2 ( 1 . 6 ) 
T=T ( z , N s , P , 5 ,BC ) 

s A 

where P g is pressure in the reaction zone and BC denotes boundary 
conditions (1.4). The solution (1.6) is considered in the flamelet 
approach as an instantaneous relations between the reactive 
species mass fractions and temperature from one hand and mixture 
fraction, scalar dissipation at the stoichiometric surface, local 
flow velocity and pressure from the other. 

Additional simplifications are possible for the low-Mach 
number combustible flows. The role of U 2 /2 term in (1.6) can be 
neglected. The pressure in the flamelet model equations can be 
treated as some constant. As the result one has: 


C =C (z, N S ,BC) 
a a ' 

T=T ( z , N S ,BC) 


(1.7) 


It is seen that flamelet model equations are splitted fully from 
the hydrodynamical ones in this case. This feature allows to 
perform the flamelet model calculations before the starting of the 
hydrodynamics one. The flamelet model eqs.(1.5) are solved for 
different values N s . The value of N s is considered as input 
parameter in these calculation. The obtained solutions (1.7) are 
collected in some database (flamelet library) in a parametric form 
on z and N s . Further to obtain the mean values of temperature, 
density and reactive species mass fractions the joint probability 
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density function (pdf) p(z,N s ) is needed. So the averaging 
procedure is reduced to the pdf model for passive scalar field 
only . 

In the current study, we tried to generalize flamelet library 
concept for the case of compressible jet flames with relatively 
small pressure gradients (it is such kind of H 2 /air diffusion 
flames [24,25] were proposed as the test cases) . That is why 
additional simplifications were adopted: i) the role of the 

pressure fluctuations on the combustion chemistry was ignored; 
ii) correlation between flow velocity and mixture fraction 
distributions was applied in eq. (1.5) in a simplified form, which 
is approximately valid for unconfined jets [28] : 

V-m - ^ d.8) 

u - u 


where lAu are the mean flow velocities of the air and fuel 
respectively, 13 is some exponent which was chosen as /3«1/Sc t (Sc t 
is the turbulent Schmidt number) . 

Such treatment allowed us to split flamelet model equations 
from the hydrodynamical ones and to apply flamelet library 
concept. The flamelet model eqs.(1.5) were solved for different 
values of N s and P . The obtained solutions were collected in 

S 

flamelet library in a parametric form on z, N s and pressure P : 

s 


C =C (z, N s , P , BC) 
a a ' ' s' 

T=T ( z , N s , P , BC) 

S 


(1.9) 


To obtain the mean values of temperature, density and reactive 
species mass fractions the joint probability density function 
(pdf) p(z,N s ) was used since the role of the pressure fluctuations 
on the combustion chemistry was ignored in the current 
calculations . 

It is convenient to demonstrate chart of the possible flamelet 
model solutions using some particular example. Such an example is 
given in Fig. 7a ( H 2 /air diffusion flame, conditions are given in 
Fig. 3). Here water mass fraction distributions are plotted vs z 
for different values of scalar dissipation N s . 

It is seen that flamelet model equations give the chemically 
equilibrium solution in the case N s =0 sec -1 (see also eq. (1.2) : 
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R a (C lf . . . ,C j( P,T) =0 a-1 , . . . , J 

Grow of scalar dissipation N s increases nonequilibriumness of 
the chemical processes (due to the increasing of the fluxes of 
reagents into the reaction zones) . Such developing 

nonequilibriumness of the water mass fraction distributions is 
seen in Fig. 7a (distributions for N s =10, 500, 969 sec” 1 

respectively) . 

When value of N s becomes too high (for the considering here 
example N s ^970 sec -1 was found) the flame extincts and the flamelet 
model gave mixing solution: 

C =(C F -C A )z +C A a=l , . . . , J (I. 10) 

a a a a 


«=1, . . . , J 


accompanied by negligibly slow oxidation. This solution is 
obtained from the flamelet model equations when the chemical 
source terms are kept to zero in (1.2) . 

The scalar dissipation value at flame extinction is referred 
to as critical value of scalar dissipation N . The burning 
solution exists only if N s < N . This means that the flame 

cr 

extincts when the mixing rate becomes too high compare to the 

fuel and air consumption inside the reaction zone due to the 

. . . . . . *) 
limitations of the finite chemistry . The calculations 

demonstrate that N cr is a pure chemical characteristic depended 

only on the fuel detailed oxidation chemistry and boundary 

conditions (1.4): 

N = N (P , Ta, Tf , kind of fuel) 

cr cr s 

where Ta and Tf are incoming air and fuel temperatures. 

The typical values of N for diffusion combustion of 

cr 

different fuels in air at room conditions (P=0.1MPa, T=300K) are 

summarized in Table 1. 

Table 1. 


Fuel 

H2 

C 3 H8 

CH4 

Ncr, sec' 1 

121 

51 

17 


’’Features of the FL predictions for the flame ignition/extinction 
regimes are discussed in Sec. III. 2 in more detailes . 
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Value of N increases with the increasing both incoming 

cr 

temperatures and pressure [22] . 

Due to the discussed here switch on/off property of the 
flamelet model solutions the range of the scalar dissipation 
variation in the flamelet model calculations is restricted by the 
range [0,Ncr) . For the higher values of N s the pure mixing 
solution (I. 10) can be used to calculate the mixture composition 
and thermodynamics properties such as density and enthalpy. 


1.2 AVERAGING PROCEDURE AND PHYSICAL INTERPRETATION 

It is seen that flamelet model predicts dependence of the 

reactive scalars on only two characteristics of the scalar field 
i.e. mixture fraction z and scalar dissipation at the flame front 
N s . The joint pdf for reactive species mass fractions C , 

temperature T, z and N s has the form: 

j 

p(z,N s ,C, . . . ,C jf T) = p(z,N s ) 8 ( T - T f 1 ) IlSfC-C^ 1 ) 

a = i 

where T fl =T(z,N s ) and C ri =C (z,N s ) are the solutions of the 

a a 

flamelet model equations (1.5). So, only p(z,N s ) requires 
additional modeling. 

To approximate the joint pdf of mixture fraction and scalar 
dissipation and to obtain the averaged values of the reactive 
species mass fractions the pdf approach [1,17] is used. Its 

features are as follows. 

The scalar field is considered as to be divided into two 

intermittent parts: i) turbulent mixing layer (0<z<l) ; ii) flow 

outside the turbulent mixing layer (z=0) . The role of the pure 
fuel flow ( z =1 ) is neglected. 

It is expected that scalar dissipation at the flame front N s 
and mixture fraction z are statistically independent inside the 
turbulent mixing layer. The role of the scalar dissipation 
fluctuations is neglected and its value conditionally averaged 
over the time moments when the turbulent mixing layer is observed 
in a given point is used. 

The Favre joint pdf of mixture fraction z and scalar 
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dissipation at the stoichiometric surface N s is considered in a 
following form: 


p/p- p(z,N s )= (1-?) 5 (z) 5 (N s ) + ■y-p t (z) 6 (N s -N®) 


(I .11) 


where 7 is the intermittency factor; p is density; p t is the 
mixture fraction probability density function in a turbulent 
mixing layer; 5 is the Dirac function. The intermittency factor j 
is calculated using approximate relation [1] : 
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'1.31/ (1+cr 2 / (z) 2 ) 

1 


if 

if 


cr/z >0.555; 
cr/z <0.555; 


( 1 . 12 ) 


where z=pz/p is Favre averaged mixture fraction and cr 2 = pz " z " /p is 
mixture fraction variance. The approximate relation (1.12) is 
based on the assumption that the fluctuation intensity inside the 
turbulent mixing layer (a^/z^) is some fixed number. Its value was 
obtained from the consideration of eigenvalue problem for the pdf 
equation (cr t /z «0 . 555) . This approximation was verified in 
[1,29,30] using various experimental data. 

There were the following reasoning to neglect scalar 
dissipation fluctuations in averaging procedure. The solutions of 
the flamelet model problem (both obtained analytically and 
numerically) predicts relatively weak dependence of the reactive 
species mass fractions on scalar dissipation at the stoichiometric 
surface N s . For example, the OH mass fraction depends on N s as 

s 1/3 

(N ) [1], [13]. That is why, the using of this simplification 

does not lead to the significant errors in averaged distributions 
of the reactive species. 

Obtained in [1] self-similar solutions of pdf equation are 
used to approximate the mixture fraction pdf P t (z) in the 
turbulent mixing layer (0<z<l). It is adopted that: 


i) . p t has the gaussian form in the non- intermittent 
the mixing layer (deeply inside the mixing layer) : 

P t = 1 exp 

y 2n a 


( z - z ) 
2 a 2 


{7=1) part of 
(1.13) 
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ii) . p t has the form of Airy function (Fig. 2) in the intermittent 
(y< 1) regions: 

p t = 1 '^° 4 Ai(l. 788 2.338) (1.14) 

t t 

where z^= z/y is the mixture fraction value conditionally averaged 
over the moments when the turbulent mixing layer is observed in a 
given point . 

The conditionally averaged value of scalar dissipation at the 
stoichiometric surface is approximated as: 



N z = zs 


H Z = Zs 


(1.15) 


where N|z = zs and 3 r|z = zs are the mean value of the scalar 

dissipation N and intermittency factor k calculated under the 
condition that mean value of mixture fraction z=zs. The 

conventional for the turbulence modeling approximation for the 
mean scalar dissipation N is used: 

Y 2 

N=0 . 07 (1.16) 

where K is turbulence kinetic energy, vt is eddy viscosity. In 
such a treatment only turbulence kinetic energy K, eddy viscosity 
vt, mean mixture fraction z=pz/p and its variance c r 2 = pz " z " /p are 
needed to calculate the pdf in any given point of the flowfield. 
These turbulent mixing characteristics are calculated using 
conventional semi-empirical transport equations of the turbulence 
modeling . 

The formulated here flamelet approach together with its 
averaging procedure has the following physical interpretation. The 
typical values of the mixture fraction stoichiometric values z 

s 

are small enough (z«0.03 for H /air flames; z -0.05-0. 06 for 

s 2 s 

different hydrocarbons /air flames) . This means that the flame 
front is located close to the outer boundary of the mixing layer. 
In this region, turbulent large-scale movement governs the mixture 
fraction large-scale fluctuations associated with the 
intermittency phenomenon. Small thickness of the reaction zone 
allows to consider it as to be "frozen" into this large-scale 
turbulent movement and to dropped out large-scale turbulence 
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influence on the chemical processes inside reaction zone. The role 
of the large-scale fluctuations is taken into account only through 
the pdf in the averaging procedure. At the same time the 
small-scale turbulence can influence the chemical process in the 
reaction zone. This influence is taken into account through the 
parameter N^, since N s characterizes fluxes of substances to the 
reaction zone [2] . The flamelet model is based on the assumption 
that this is the main cause which is responsible for the influence 
of turbulent mixing on the local nonequilibriumness of combustion 
chemistry . 


1.3 COUPLING WITH COMPRESSIBLE CFD SOLVERS 


The following procedure was proposed for the flamelet library 
incorporation into the compressible flow hydrodynamics solver. The 
"effective" heat capacities CP ft of the reactive species are 
introduced in the same manner as it was done in [31] : 

.T 


CP a = fcp a dT/ (T-To) 

T o J 


(1.17) 


Using (1.17), the total mixture enthalpy can be written as: 

j -> -> 

H = CP- (T-To) + V C Ah + (U-U)/2 

a a 

0C= 1 


;i . is; 


j 

where CP=^ CP a C a is the "effective" heat capacity of the mixture. 


a= i 


Let us introduce two additional "effective" parameters. 
"Effective" heat capacities ratio T which is defined as: 

r=l/(l- R/(CP-p) ) 


(1.19) 


and "effective" heat of mixture formation defined as: 

j 

0=^ C a Ah a -CP-To (1.20) 

a= i 

where 4 is the mixture molar weight and R=8.31 J/ (mol K) is 
universal gas constant 

Using the thermal equation of state for mixture P=pRT/u and 
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eqs. (1 .19) , (1 .20) one can obtain from (1.18) the following 
"effective" form for the total enthalpy: 

-> -» 

H=r/(r-l) P/p + Q +(U-U)/2 ( 1 . 21 ) 

Multiplying eq. (1.21) by density and Favre- averaging one has: 

p H = [r/(r-l)] P + pQ + pU 2 / 2 + pK (1.22) 

where K is turbulent kinetic energy. 

Additional simplifications were adopted i.e. the correlation 

(P'T') and K in (1.22) were neglected. As the result one has: 

p H = [r/ (r-1) ] P + p Q + pU 2 / 2 (1.23) 

It is seen from (1.19), (1.20) that values of r/(r-l) and Q 
depend only on reactive species mass fractions and temperature. 
That is why, they can be obtained from the flamelet calculations. 
To obtain the mean values of r/ (T-l) and Q the averaging procedure 
of Sec . 1 . 2 can be used. The value of H is obtained from the 
Favre -averaged energy conservation equation in its conventional 
form. As the result, eq.(1.23) give the relation between mean 
values of density p, pressure P and flow velocity U where only two 
parameters (T and Q) are needed from the flamelet model 
calculations. Of course such a simplified procedure does not allow 
to calculate all the mixture thermodynamics properties (for 
example, the local speed of sound) however it allows significantly 
reduce computational expenses. If mean mixture composition is 
needed in some cross sections the whole flamelet library should be 
used. 
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Part II. CFD TESTS FOR SUPERSONIC JET FLAMES 


CFD tests of the flamelet model capabilities were done for two 
test cases [24,25] where H 2 jet combustion in supersonic air flow 
was studied experimentally. The goals of the tests were: 

- to estimate computational expenses for FL realization in CFD; 

- to examine accuracy of model predictions. 


The flamelet library concept of Sec.I.l was applied (FL 
approach) . At the first step, flamelet equations (1.5) were solved 
parametrically and obtained solutions were collected into the 
flamelet library. At the second step, flamelet library was used 
together with the appropriate CFD solvers for the flowfield 
calculations using procedure of Sec. I. 3. 

Additional series of calculations were done using 
"quasilaminar " combustion model together with the same CFD solvers 
(QL approach) to obtain the reference point for comparison of FL 
approach computational and physical capabilities. The mass, 
momentum and energy conservation equations, in QL approach, were 
solved together with the averaged conservation equations for the 
reactive species where the role of the reactive scalars turbulent 
fluctuations was neglected: 


|-pU C = p 
dxT i a 5x K 

i i 


(-§^ + D )H“ + P R a <f ’P’ e , <V 

V t ' i 


(II. 1) 


a=l , . . . , J 


where are the Cartesian coordinates; the chemical source terms 
are postulated in Arrenius form for the mean values of species 
mass fractions 2 , temperature T and density p. 

Both FL and QL series of calculations were done using the same 
initial and boundary conditions, model of turbulence and 
approximations for the detailed kinetics and thermodynamics 
properties . 

The detailed hydrogen oxidation chemistry was approximated by 
the Miller-Bowman kinetics scheme [32] . The thermal NO formation 
mechanism was taken into account also. The resulting detailed 
kinetics model included 21 reactions between 11 species (H 2 , 0 2 , 

H 2 0, H, 0, OH, H 2 0 2 , H0 2 , N 2 , N, NO) . It is given in Appendix C. 

The species specific enthalpies h a = fc Po( dT + Ah a (To) were used in 
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7 

a form of polynomial approximations h = A + V A. (T/1000) 1 

OL I OL Zj i CL 

i = 1 

taken from [33] . Here Ah^ are the species heats of formation at 
reference temperature To= 298. 15K. The polynomial coefficients are 
given in Appendix C also. 

The account about details of computations and obtained results 
are given in Sec. II. 1 and Sec. II. 2. These results and their 
detailed analysis were presented in previous interim Reports [27] . 
Additional methodological tests of the computations accuracy, 
which were done at the final step of the research, did not change 
results and conclusions of [27]. So, Sec. II. 1 and Sec. II. 2 mostly 
compiled [27] . Brief summary and comparative analysis of the 
results are presented in Sec. II. 3 together with our suggestions 
concerning ways for the further model improvement . 
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II .1 BEACH COAXIAL EXPERIMENT [24] 


DESCRIPTION OF THE TEST CASE 

The sketch of the Beach et al . test case [24] is given in 
Fig. 3 together with the nozzle exit conditions (flow parameters 
and gas composition) . The hydrogen was injected through supersonic 
axisymmetric nozzle with the Mach number M =2. The hot air was 

H2 

obtained by burning of hydrogen in air, replacing the oxygen and 
expanding through supersonic nozzle with the Mach number M =1.9. 

air 

The hydrogen injector tube had external diameter d=0. 009525m with 

j 

a lip thickness 0.0015m. The air nozzle free stream diameter D was 
0.0653m. 

FLMELET LIBRARY GENERATION 

The boundary conditions for the flamelet equations were 
adjusted according to the data presented in Fig. 3. The influence 
of the pressure variation inside the flowfield on the combustion 
chemistry was neglected and flamelet model equations were 
calculated for a fixed value of pressure P =0.1MPa. The value of 

S 

exponent £ in approximation (1.8) was chosen as £ = 1.25 which 
corresponded to the value of the Sc t =0 . 8 . 

The flamelet eqs.(1.5) were solved using time-relaxation code 
FLSLV (its description is given in Appendix B) for different N s to 
cover whole range of the possible flamelet model solutions from 
N s =0 up to N s =N cr . The obtained distributions of reactive species 
C a (z,N s ) (a=l,...,J) and static temperature T(z,N s ) were 
introduced into the flamelet library. The whole flamelet library 
was calculated using exponential grid consisted of 1=81 points 
with grid points clustering near z=0 boundary. 

Computational strategy and convergence 

The calculations were started from lowest value of N s =N (1) (it 
was chosen as N (1) =0.001 sec -1 ). The chemically equilibrium 
solution for reactive species and temperature was used at this 
step as the initial approach. Further the converged solution for 

■ s (1 ) 

the first value of N =N was used as the initial approach in 
calculations performed for the next value of N s = N <2) and so on. 

The increments in time SC control was applied for choice of the 
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optimal pseudo-time step x and to fasten the convergence to steady 
state solution (see detailed in Sec.B.I of Appendix B) . The norm 
for SC was calculated at each time step as: 

e j = mgx max{mod (SC a ) (C a ) ^ _1 } 

where maxima was got through all species (a=l, ...,11) and all grid 

points z i (i=l,...,I). The time relaxation was stopped when the 

value of e j became lower than c =10~ 6 . The obtained solution was 

0 

considered as converged. Typical residual norm e J and time step r J 
behavior in course of flamelet library generation are shown in 
Fig. 4. The required number of iterations to obtain the converged 
solution was 80 at the first step of the flamelet library 
generation and then it rapidly decreased up to the value 5 . The 
required iterations number chart during the flamelet library 
generation is given in Fig. 5. 

Accuracy of the computations 

Two tests were done to estimate the accuracy of the obtained 
flamelet library. 

The first one was the extrapolation of the obtained numerical 
solution to zero-length grid step (h-»0) . For this purpose the 
additional methodological calculations for total grid point number 
1=41 and 1=161 were performed for three selected values of scalar 
dissipation N s (N s =0 . 01 ; 100; 900 sec 1 ) . The obtained reactive 
species profiles were compared with those which were obtained at 
"basis" grid 1=81. To avoid errors associated with the application 
of extrapolation procedure to the case where the numerical grid 
was nonuniform the fine grid cells were generates by dividing of 
the rough grid cells in half strictly. The results of calculations 
obtained for three grids ("basis" 1=81; "fine" 1=161 and, for the 
control, "rough" 1=41 ) were used in extrapolation of the solution 

to the zero- length step solution C (Rch ! The following norm for the 
accuracy of solution was introduced: 

ERROR <Rch) = max max mod [ (C ( 81) -C <Rch) ) /C <Rch) ] 
a .,s ai oa on 

N i 

where maxima was got through three control solutions corresponded 
to N s =0.01, 100, 900 sec 1 and all grid points z. of the "basis" 
grid (1=81) . The reactive species mass fraction values lower than 
10 C a were ignored in this estimation. It was found that 
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ERROR (Rch) was lower 0.3% for "effective" heat capacities ratio r, 
a 

"effective" heat of formation Q and all substances except H^. 

For HO its value was 1%. However this feature can not influence 

2 2 

the accuracy of the flamelet library as a whole since the HO 
mass fraction was too small ( < 1 0 -6 ) and could not influence total 
mixture thermodynamics properties. The examples of the 
extrapolation procedure for maximum concentrations values of two 
substances (H, OH were selected) are given in Fig. 6. They are 
plotted vs l/I where I is the total number of grid points. The 
solid line in Fig. 6 corresponds to the mean root square linear 
approximation . 

The second test of the flamelet library accuracy was 
associated with the fact that the CFD will require to obtain the 
reactive species profiles at interim values of N s . To investigate 
the accuracy of the interpolation the additional calculations were 
performed for interim values of N° +i/2 = (N% N s +i )/2; (r=l, ...,37) 

where N% N* +i are scalar dissipation values corresponded to the 
particular solutions which were included to the flamelet library . 
The results were compared with those obtained by the linear 
interpolation of the flamelet library data between solutions at N s 

r 

and at . The accuracy of the interpolation was defined as: 


ERROR ( 1 nt) = m ax 

(X s 

N 

r + 1 / 2 


max mod [ ( C 


(N) ( int) . /r (N) , 

ai on ' on J 


where maxima was got through all interim solutions 


points z of the "basis" grid 


, ( N) 


flamelet calculations at N S =N S 


(1 = 81); C'‘” denotes 

and C ( lnt) denotes 

r+l/2 ai 


and all grid 
results of 
results of 

interpolation between solutions containing in the library. It was 

lower than 0.5% for r,Q and species 
10" 6 . The ERROR^ int) value was found to 

concentration of 
influence total 


found that ERR0R^, lnt) was 


concentrations higher than 
be about 10% for radicals 
these radicals was small 
mixture properties. 


H 0 and HO but the 

2 2 2 

(<10 ) and can not 


Obtained results and computational expenses for flamelet library 
generation. 

Total flamelet library included 38 particular solutions in the 
range of N s variation from 0.001 sec -1 up to N = 970 sec" 1 . 

cr 

Additional time was required (8 particular solutions) to adjust 
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the value of N with the accuracy 0.1%. Pure mixing solution 
(1. 10) was postulated for N s iN_ r . Examples of the obtained 
distributions in parametric form on N s are given in Fig.7a-c for 
stable species (H 2 0, H 2 , 0 2 ), i n Figs.8a-c for main radicals 

(OH, H, 0) and in Figs. 9a, b for "effective" parameters r,Q. 

The flamelet library generation was done using conventional PC 
AT 486DX2 /66MHz computer. The computational expenses required for 
the generation of the flamelet library are summarized in Table 2 . 

Table 2 


Number of species 

11 

Grid in z direction 

81 points 

CPU time per iteration per grid point 

0 . 0092sec 

Total time for flamelet library generation 

607sec 

Total number of calculated 


particular solutions 

46 

Number of particular solutions 


included into the library 

38 

Required memory for library storing 

0 . 8Mb 


Tests of results sensitivity 

The methodological tests were done to investigate the 
sensitivity of the flamelet model calculations to the choice of 
the detailed kinetics scheme for conditions of the Beach test 
case. The solution of the flamelet model equations for three 
selected values of scalar dissipation N s (N s =l, 100 , 800 sec” 1 ) was 
obtained using detailed kinetics for H 2 oxidation proposed by 
Warnatz in [34] (it is presented in Appendix C) and compared with 
the results obtained by Miller-Bowman scheme. The examples of 
obtained results are given in Fig. 10 for two substances (HO and 
OH) . We have not found any significant influence of the detailed 
kinetics approximation on the results of calculations for main 
reactive species. For example, the difference for H 2 0 
concentrations was about 2% and for maximum OH concentrations it 
was about 15%. 

The additional methodological test was done to investigate the 
sensitivity of the results to the exponent /3 in the adopted 
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correlation (1.8). The parameter £ was varied in the range 

£=0.5 -2.0. The water and hydroxil radical mass fractions 

(C H 0 > C 0H ) distributions obtained at scalar dissipation value 
2 

N s =50 sec and different values of (3 are presented in Figs. 11a, b. 
It was found that the sensitivity of the results to the £ 
variation is relatively small (less then 1% for main stable 
species and less than 12% for radicals ) . It is compatible with 
the sensitivity to the adopted detailed chemistry approximation. 

FLOW FIELD CALCULATIONS 

Adopted simplifications and system of equations 

The assumption about H 2 jet in co-flowing infinite air stream 
for flow hydrodynamics was adopted. The role of the air flow 
mixing with the ambient air was neglected. The features of the 
flowfield in the vicinity of the hydrogen nozzle exit lip were 
taken into account only through the initial conditions in the 
initial cross section of the computational domain in cross-section 
x/d.=0 . 33 downstream the injector. It was expected that the 
flowfield can be described by the parabolized approximation (PNS) 
of the 2-D Favre- averaged conservation equations. 

A special approach concerning governing system of equations 
was adopted in current study to provide the stable marching 
calculations of PNS equations in slightly subsonic regions which 
can arise inside the mixing layer for the conditions of Beach test 
case [24] . For this purpose the procedure of PNS equations 
regularization proposed in [35-37] was applied. The term with 
longitudinal pressure gradient in x-momentum equation was 
multiplied by the parameter u, which was expected to be function 
of the local Mach number estimated on the longitudinal velocity 
component U. The Cauchy problem with initial data for PNS 
equations in subsonic regions (M <1) is well posed at the 

condition [37] : 

..2 

u < M 

X 

In the pure supersonic regions (M >1) parameter w was equal to 1. 
The influence of the rejected part of the longitudinal pressure 
gradient (1-w) on the solution was neglected since the 

pressure gradients are small for the conditions of Beach test 
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case. So only downstream marching without global pressure 

. . * ) 
iterations was applied. 

The system of the regularized PNS equations had the form: 
continuity equation: 


d (pU) 


d (pV ) 
ay 


pV 


3x ay y 

X-momentum equation: 


= 0 


a (pUU) 


a (pw) 

3y 


pUV 


dx ay y 

Y-momentum equation: 


a_ 

ay 


- au 

^tfly 


pvt au 
y 3y 


- u> 


ap_ 

dx 


a (puv) 

dx 


d(P W+p) | pW 
ay y 


a f 2 - 
= ay la 


av _ v]] + 2pvt rav _ v 1 
a y ~ yJJ y l a y y ; 


Energy conservation equation: 
apufi ^ apVH ^ pvfi _ a (pvt aH 


dx 


3y 


y 


ay l Pr t 5 y 

a - 


l pi^t afi 

y P^ t ay 


ayP v t 


[u 2 ^[ 2 av _ 2 ]] 
l ay 3 ( ay yj J 


(II. 2) 


(ii. 3 : 


(II .4) 


(II. 5) 


Here x,y are the longitudinal and transverse axe of the 
coordinate system; U,V are the components of the velocity vector; 
p is density; P is pressure; K is static enthalpy of the mixture 
taking into account formation enthalpies of the mixture 
components; H is the total enthalpy which is defined as 
H=K+ (U 2 +V 2 ) /2 ; v t is the eddy viscosity; Pr t is the turbulent 
Prandtl number (Pr t =0.8); the upper symbols (~) and (”) denote 
Favre and time averaging respectively; o> is regularization 
factor . 

The Secundov's one-equation turbulence model "i^- 90" was used 
to calculate the eddy viscosity i> t . Here its general form was 
reduced to the following parabolic equation: 


dpUy t dpVvt pVv t 

dx + dy + y 


d_ 
dy | 


p(c Vt + v) „ 

K i dvt 


( P< c , l ’ t * •'» -fv - ) 

ay + - c.P*” 1 


(II. 6) 


+c 


- 2 G 


*)The validity of the marching PNS calculations for the test case 
[24] was checked at the final step of the study. The full NS solver 
FNAS2D (see Sec.B.3 of Appendix B) was used for calculations of 
the Beach test case. The minor difference in results was found. 
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where G= 


; a is local speed of sound; c t = 2; c 3 =0.7; c 4 =5; 


c = 0.2 

2 (30 vt) 2 + y 4 G 2 


; v is kinematic viscosity. 


To obtain values of the mean z and variance cr‘ 


mixture 


fraction values and mean scalar dissipation N the following 
equations were solved: 

mixture fraction z=pz/p transport equation 


9_puz + 3_pv£ + P Vz = 2_f _PELdz ) _i_ pv _ t az. 

dx P 9yP y Sy[ Sc t 3y J y Sc t dy' 


mixture fraction variance <x = pz"z"/p transport equation 


dpUcr 2 apVo- 2 pVcr 2 
3x + 3y + y 


+ 2 -^ 
Sc 


a T pi^t ao- 2 ] i pt^t da‘ 
ay [sc t ay J + y sc t ay 

Vt ( dz ) 2 a pKcr 2 

Sc t [ 3y J P i vt ; 


the turbulent kinetic energy K balance equation 


apUK apVK pVK _ d 


dx + 5y + y 


= 7irr\PK vt 7F7 + TT PK vt 7577 + P vt 


where k 2 = 1.4; /3 i = 0.14; /3 2 =0.1; Sc t =0.8. 

The values of z, cr 2 were used to calculate the mixture 
fraction pdf p(z) using formula (1 . 11) - (1 . 14) of Sec 1.2. 

The conditionally averaged value of scalar dissipation at the 
flame front if* was calculated by approximation (1.15) using the 

v 2 

mean scalar dissipation N=0.07 — and intermittency factor y 
distributions . 

The particular solution of the flamelet model problem at N S =N^ 
was obtained from the flamelet library by the linear interpolation 
between neighboring solutions for N (k) and N (k+1) , where 

-_(k) - T s -_<k+l) .. . . _ 

N and linear interpolation on z. The obtained 

solution for the "effective" heat capacities ratio F and effective 
heat of "formation" Q were averaged using calculated p(z) . The 
averaged values of Y] (T-l) and Q were used for the closure of the 
governing system of conservation equations ( II . 2 ) - (II . 6 ) using 
procedure outlined in Sec 1.3. 
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Boundary and initial conditions 

The computations were done for the rectangular domain 
presented in Fig. 12. The calculations were started in the cross 
section x/d^O.33 and they were stopped in the cross section 
x/d j =30. The upper boundary of computational domain was at y/d.=2 
position. The no-reflection conditions were posed at the upper 
boundary of the computational domain. The symmetry conditions were 
posed at the axis of symmetry. 

The parameters profiles adopted as the initial conditions at 
cross section x/d^O.33 are presented in Fig. 13. The initial 
profiles for the longitudinal velocity and turbulent 
characteristics were chosen based on estimations of the boundary 
layers thickness on the external and internal sides of hydrogen 
injector and additional turbulence production in the wake 
downstream nozzle lip. The non-dimensional initial distributions 
of the longitudinal component of the flow velocity U (o) =U/U ; 

/ \ H2 

eddy viscosity i>° =v / (U d) ; turbulent kinetic energy K ( °=K/U 2 
as well as initial distributions of z, cr are given in Fig. 13. The 
transverse component of the velocity V was expected to be zero. 
The species mass fractions (H 2 , C> 2 , ly}, N 2 ) were expected to be 
constants in the inner and outer flows and their values were 
chosen in accordance with data of Beach [24] . In the intermediate 
region they were postulated as linear functions of mixture 
fraction z. All other species concentrations were equaled to zero. 

The total enthalpy distribution was expected to be uniform for 
pure H 2 and pure air flows. In the intermediate region it was 
expected to be linear dependent on z. 

Minimum information concerning nozzle configuration and 

initial distributions was given in the original paper of Beach et 

al . [24] . However we have tried to make an indirect estimation of 

adopted distributions for U, H, and z validity. For this purpose 

we calculated the pitot pressure distribution in the initial cross 

section using the adopted initial distributions and compared it 

with the experimental data of Beach (Fig. 14a) . The validity of the 

adopted profile for eddy viscosity was approximately estimated 

from the correspondence of the centerline H mass fraction 

2 

distribution obtained in FL approach calculations with the 
experimentally measured by Beach (Fig. 14b) . It is seen that basis 
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features of the experimentally measured distributions were 
reproduced correctly. 

CFD solver and numerical grid 

The regularized PNS equations were solved using modified 
version of the marching code SUPNEF. The code SUPNEF is based on 
explicit finite difference method which is the generalization of 
the well-known steady analogy of Godunov method [38] for the 
steady supersonic flows. The code description is given in Sec.B.2 
of Appendix B . 

Code modifications were connected with the incorporation of 
the regularization procedure for slightly subsonic regions in 
accordance with [39] . For this purpose the characteristic 
relations for the inviscid part of the regularized PNS equations 
were used for the solution of two flows interaction problem based 
on the assumption that intensity of main discontinuities is small . 
The obtained solutions were used to approximate the convective 
fluxes on the cell boundaries in subsonic regions (regularization 
factor (j<l) . In supersonic region (w=l) these relations coincide 
with usual relations for steady supersonic flows interaction 
problem. The following relation for the regularization factor u 
was used in real calculations: 



The calculations were done using adaptive grid. The grid 
adaptation was realized in accordance with spring analogy [40] . In 
each cross-section all grid nodes were supposed to be connected by 
springs with the stiffness proportional to the gradient of Mach 
number. The nodes positions were determined in accordance with 
springs system equilibrium conditions. 

Methodological tests and computational expenses 

The methodological calculations were performed using 50 , 100 
and 200 computational cells in cross sections for both FL and QL 
approaches. The role of the grid nodes number variation on the 
calculated distributions of reactive species mass fractions and 
turbulent mixing characteristics is given for both FL and QL 
approaches in Fig. 15a, b respectively for x/d^S.26 cross section. 
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All the final FL and QL computations were performed using 
adaptive grids containing 100 computational cells in each cross 
section. They are given for both approaches in Fig. 12. 

The computations were performed using workstation HP 9000/735. 
The CPU time requirements are given in Table 3. 


Table 3. Computational expences for flowfield calculations 



QL approach 

(21 reactions, 
11 species) 

FL approach 

(21 reactions, 
11 species) 

Number of grid cells 

100 

100 

CPU time per grid cell 

0 . 00162sec 

0 . 00097sec* ) 

Total CPU time of computation 

424sec 

260sec 


*) CPU requirements connected with fiamelet library interpolation, 
calculations of p(z) and averaging were 0. 00017sec/cel 1 


RESULTS OF COMPUTATIONS 

The obtained H 2 0 mass fraction contours are given in Fig. 16 
for both FL and QL approaches. It is seen that fiamelet model 
predicted transition from mixing to burning regime in the vicinity 
of x/d j= 3 cross section. Here the conditionally averaged value of 
the scalar dissipation at the flame front isr* became lower its 
critical value N (Fig. 17). 

cr 

The obtained Mach number contours are given in Fig. 18 for both 
FL and QL approaches . It is seen that mixture ignition in the 
fiamelet model calculations is accompanied by the sharp increasing 
of the released heat and generation of the weak compression wave 
and slightly subsonic region. It is seen that QL approach predicts 
more smooth heat release increasing in the mixing/burning 
transition region. 

The cause of such difference between FL and QL predictions is 
seen from the consideration of Fig. 19, where the C> 2 mass fraction 
profiles obtained slightly upstream and downstream ignition point 
are presented for both FL and QL approaches . Fiamelet model 
generates solution which corresponds to the approach where all the 
fuel, which initially penetrated fuel-lean part of the mixing 
layer, reacts with the oxygen in narrow vicinity of the ignition 
point . The QL approach generates solution where fuel and oxygen 
consumption is much more weak. It is observed only in the regions 
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with near- stoichiometric mixture composition, which is more 
realistic in nature. Fortunately, the ignition distance for the 
conditions of Beach test case was small enough and the amount of 
the H 2 penetrated the fuel -lean part of the mixing layer was about 
1.6% of its total mass flow rate. That is why the disturbance of 
the flow field was small also. Such flamelet model behavior is not 
important for the flames with the short ignition delay length or 
stabilized in the vicinity of the fuel injectors. However it can 
be serious disadvantage of the model in the case of flames having 
large ignition delay. 

The obtained in FL approach contours of the turbulent mixing 
characteristics (mean and variance mixture fraction, turbulent 
kinetic energy) are given in Fig. 20. 

The obtained in FL approach distributions of the averaged mass 
fractions of H 2 0, H 2 , C> 2 and N 2 are given in Figs.21a-d by solid 

lines together with the data of Beach for four cross sections 
(x/dj= 8.26,15.5, 21.7 and 27.9). Good correlation between FL 

predictions and experimental data is seen. As a rule, FL 
predictions of species distributions were possible with accuracy 
better than 20%. Much discrepancy («25%) was observed only in 

predictions of H 2 0 peak value and for only one test section 
(x/dj = 8 . 26) . Such discrepancy can be explained taking into 
account that the estimated error of measurements was higher than 
15% (due to mentioned in [24] possible mixture reacting inside the 
sampling probe) . It can be attributed also to the kind of averaged 
values (Reynolds or Favre) measured by sampling technique. 

Results of the QL predictions are given in Fig.21a-d by dotted 
lines. It is seen that the FL predictions are closer to the 
experimental data as a whole compare to the QL ones . The QL 

approach gave significant overprediction of 0 2 mass fraction in 
fuel-rich regions and displacement of the H 2 0 peak location in 

comparison with experimental data. 

It was concluded that FL approach gave satisfactory result for 
Beach test case [24] . 
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I I. 2 BURROWS -KURKOV EXPERIMENT [25] 


DESCRIPTION OF THE TEST CASE 

Fig. 22 shows scheme of setup and test conditions of 

Burrows -Kurkov experiment [25] . The test section was rectangular 
duct having the constant width (0.051m). The air supply duct had 
0.089m height. Hydrogen was injected parallel to the vitiated air 
flow. It was injected with a sonic speed through the 
two-dimensional slot located at the backward step in the initial 
cross section. Slot height was h=0.004m. Lip thickness at the top 
of the step was 0.76*10 _3 m. Test section total height expanded 
linearly from 0.0938m in the initial cross section to 0.105m at 
the exit cross section. Composition measurements were done at the 
exit plane of the test section located at x=0.356m downstream the 
injector location. 

FLAMELET LIBRARY GENERATION 

The procedure of the flamelet library generation was basically 
the same as that used for the Beach test case. The boundary 
conditions for the flamelet equations were adjusted according to 
the data presented in Fig. 22. The thermochemistry approximation of 
Appendix C was used for species enthalpies and detailed chemistry 
model (21 reactions between J=ll species) . Only one additional 
feature was taken into account. It was the influence of the 
pressure variation inside the flowfield on the combustion 
chemistry. The range of the pressure variation was adjusted based 
on the presented in [25] static pressure distributions along the 
duct wall («0 . 08-0 . 12MPa) . The flamelet model calculations were 
done for the pressure values in the reaction zone 
P=P s =0 . 08 , 0 . 1 , 0 . 12 MPa. The obtained solutions were united into 
the total flamelet library. 

The value of exponent 0 in approximation (1.8) was chosen as 
0=1 which corresponded to the value of the Sc =1. 

The whole flamelet library was calculated using the same as 
for Beach test case exponential grid consisted of 1=81 points. 

Total flamelet library included 171 particular solutions in 
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the range of N s variation from 0.001 sec Up tO N = N (P ) 

cr cr s 

where the later quantity was calculated for three values of 
P s =0.08, 0.1 and 0.12 MPa. The N cr (P s ) dependence was approximated 

as N_^=777 .43* (P^/0 . IMPa) -135 .419 sec 1 with the accuracy 0.4%. 

Examples of the obtained distributions of temperature T, H 2 0, 
OH radical, and "effective" parameters Q and r are given in 
Fig.23a-e in parametric form on N s . The influence of the pressure 
variation on the results of flamelet model calculations is 
illustrated by Fig.24a-c. 

The computational expenses required for the generation of the 
flamelet library at PC AT 486DX2/66MHz are summarized in Table 4. 


Table 4. 


Number of species 

11 

Grid in z direction 

81 points 

CPU time per iteration per grid point 

0 . 0092sec 

Total time for flamelet library generation 

2368sec 

Total number of calculated 


particular solutions 

185 

Number of particular solutions 


included into the library 

171 

Required memory for library storing 

3.21Mb 


The same, as for the Beach test case, norms of the flamelet 
library accuracy were estimated: 

i) relative error ERROR (Rch) associated with the Richardson 
extrapolation of the obtained numerical solutions to zero-length 
grid step; 

ii) relative error ERROR ( 1 nt) associated with the accuracy of 
the interpolation for the interim values of scalar dissipation N s 
and pressure P g using solutions contained in the flamelet library. 

It was found that ERR0R^ Hch) < 0.3% and ERR0R^ int) < 0.5% for 
"effective" heat capacities ratio T, "effective" heat of formation 
Q and for all the reactive species mass fractions having the 
maximum values higher than 10” 6 . 

Tests of results sensitivity to the detailed chemistry model 
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and adopted U=U(z) correlation were done. 

The Miller-Bowman detailed chemistry approximation was 
substituted by the Warnatz one [34] . The calculations of the 
flamelet model were repeated with the new detailed chemistry 
approximation for three different values of scalar dissipation N s 
(N s =l,30,250 sec x ) and were compared with the results obtained by 
Miller-Bowman scheme. The examples of obtained results are given 
in Fig. 2 5 for two substances ( H 2 0 and OH) . We have found minor 
influence of the detailed kinetics approximation on the results of 
calculations for main reactive species and "effective" parameters 
Q and T. 

The parameter £ in U=U(z) correlation (1.8) was varied in the 
range £=0.5-2. 5. The temperature, water and effective parameters Q 
and T distributions obtained for scalar dissipation value N s =30 
sec" 1 and different values of £ are presented in Figs.26a-d. It 
was found that the sensitivity of the results to the £ variation 
is relatively small and compatible with the sensitivity of the 
results to the adopted detailed chemistry approximation (about 1% 
for main stable species, about 2% for temperature and less than 
20% for radicals) . 

FLOW FIELD CALCULATIONS^ 

Adopted simplifications 

The flow field was expected 2-D. The role of the boundary 
layer at the upper duct wall was neglected since it was difficult 
to have satisfactory resolution for the boundary layers on the 
both walls of the duct due to limitations in operational memory 
of the available HP workstation (64Mb) . However the adopted 
approximation seems to be justified because the height of the test 
section was much greater than the slot for hydrogen injection. It 
is possible to assume that the boundary layer on the upper wall 
does not perturb the mixing and boundary layers near the lower 
wall . 

The simplest molecular transport model was applied i.e. the 
mixture molecular viscosity and diffusivity were estimated based 
on H 2 molecular diffusivity and fixed laminar Prandtl and Schmidt 
numbers Pr=Sc=0.72. 
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System of equations 

The following 2-D approximation of steady, averaged Navier-Stokes 
(NS) equations was used for the flowfield calculations: 


5F dG 
+ 

dx dy 


(II .7) 


where vectors of fluxes F,G were as follows: 
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Here: x,y are the longitudinal and transverse coordinates; U,V are 
the components of the velocity vector; p is density; P is 
pressure; H is the total enthalpy. 

The stress terms were defined as: 
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where v is kinematic viscosity; v is the eddy viscosity; 

Pr is turbulent Prandtl number (Pr =1) . 
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The Secundov' s one-equation turbulence model "v - 90" was used 
to calculate the eddy viscosity 
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Here a is local speed of sound; S is minimum distance from the 
wall . 

To obtain values of mean mixture fraction z and its variance 

2 . ~ 

<r values and mean scalar dissipation N the following transport 

equations were solved: 


apuz apuz a v t v .azl a (- . v t v . ai] 

ax + ay ax [ p Sc t + Sc ; axj + ay( p ' Sc + Sc 'ayj ; 


(II. 9) 
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where k = 1.4; /3 = 0.14; 13 =0.1; Sc =1. 

2 1 2 t 

~ 2 

The values of z, a were used to calculate the mixture 
fraction pdf p(z) using formula ( 1 . 11) - ( 1 . 14) of Sec 1.2. 

The conditionally averaged value of scalar dissipation at the 
flame front was calculated by approximation (1.15) using the 

~ K O' 2 

mean scalar dissipation N=0.07 — and intermittency factor j 
distributions . 

The following procedure was adopted to adjust the value of the 
pressure P for the selection of instantaneous flamelet model 

s 

solution. The location of the closest to the mean stoichiometric 
surface z = z s computational cell in each cross section of the duct 
was emphasized. The value of mean static pressure P in this cell 
was used as the reference pressure Ps inside the reaction zones 
for the cross section. 
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The particular solution of the flamelet model problem at N S =N* 
and Ps = P(z = zs) was obtained from the flamelet library by the 
linear interpolation between neighboring solutions for scalar 
dissipation, pressure and mixture fraction. The obtained solution 
for the "effective" heat capacities ratio r and "effective" heat 
of formation Q were averaged using calculated p(z). The averaged 
values of T/ (T-l) and Q were used for the closure of the governing 
system of conservation equations using procedure outlined in 
Sec .1.3. 

In the case of QL approach calculations the system 
(II . 7) , ( II . 8 ) was solved together with the averaged reactive 
species mass conservation equations (II. 1). 

Computational domain and boundary conditions 

The computational domain is given in Fig. 27. The left boundary 
was located in the hydrogen injection cross-section where all 
parameters distributions were supposed to be known in all opened 
parts of cross-section excluding the slot lip. No-slip velocity 
conditions were posed on the lip wall and on the lower wall of the 
duct . It was expected also that the lip temperature was fixed 
(T w =300K) . The lower wall was expected to be adiabatic with zero 
temperature gradient . The turbulent kinetic energy and turbulent 
viscosity were equal to zero at the lower wall and lip. The mean 
mixture fraction and mixture fraction variance normal derivatives 
were equal zero at the lower wall and lip. 

The upper wall was considered as inviscid with zero 
transversal velocity component. All normal derivatives, which were 
needed to estimate viscous stresses and corresponding diffusion 
fluxes on the wall, were equal to zero. 

In the exit plane of computational domain (located at x=0.356m 
cross section downstream the injector) the so-called drift 
boundary conditions with normal derivatives of all parameters 
determination from computational domain were posed. 

In the case of QL approach calculations all walls were 
supposed to be noncatalitic and species concentrations normal 
derivatives were set to be equal to zero. 

The parameters distributions at the inlet boundary were 
obtained by the following manner. The longitudinal velocity and 
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eddy viscosity distributions in the incoming air flow were derived 
from the experimental data of [25] on boundary layer thickness 
(S»0.012m). The eddy viscosity distribution inside the boundary 
layer was approximated using formula: 

yU r 

v =0.41 y U (l-exp(- — — ■) (1-y/S) exp(-y/6); 

t L <z o l/ 

where dynamic velocity U was calculated using correlation for 
local skin friction factor C for compressible turbulent boundary 
layer at flat plate: 

2U 2 /U 2 =C =0 .023 (1 + 0.7 (k - 1) M 2 / 2 ) ~°' 5 (2 / ( 1+Tw) °' 5 . 

X e f V e 

Here subscript "e" denotes parameters in core flow, Tw is the 
temperature factor, k is heat capacities ratio. 

The longitudinal velocity distribution in the boundary layer 
was obtained by integrating of the equation for the shear stress: 


Ixy 


= <v*) 


au 

ay 


where the distribution for the shear stress inside the boundary 
layer was approximated according to [41] as: 


— = U 2 (1-3 (y/<5) 2 +2 (y/S) 3 ) 

P 

The turbulence kinetic energy distribution in the boundary 
layer was calculated from the approximate "equilibrium turbulence" 
relation : 

v 

-Y |3U/ay | =0 . 3 . 

The same procedure was used for U, v and K distributions 
calculations in the exit plane of the hydrogen injection slot. The 
boundary layers on the walls were expected to be fully developed 
with the thickness equal to a half of the slot height 
(S=l/2h=0 . 002m) . The hydrogen velocity profile maximum value was 
adjusted to provide the same total H 2 mass flow rate as that 
obtained in experiment . 

The 2% velocity fluctuations were expected in free stream 
outside the boundary layers. 

The obtained non-dimensional initial distributions of the 
longitudinal component of the flow velocity U <o) =U/Ua; eddy 
viscosity ^° } =v/ (Uah) ; turbulent kinetic energy K (o iK/uI and 
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total enthalpy H (0) =H/Ua are given in Fig. 28 where Ua is the air 
flow velocity in the core flow and h is hydrogen injector slot 
height . 

The uniform step initial profile was expected for mixture 
fraction z. The transverse component of the velocity V and mixture 
fraction variance cr 2 were expected to be zero at the entry 
boundary . 

To estimate the reliability of the adopted initial 

distributions we have calculated pure mixing regime of 

Burrows -Kurkov experiment. No additional oxygen was introduced 
into the air flow in this test run and there was no combustion as 
the result. The distributions of the averaged mole fractions of 
N 2 ,H 2 and H 2 0 obtained for cross section x=0.356m are given in 
Fig. 29 together with the experimental data. The satisfactory 
agreement is seen. 

CFD solver, numerical grid, convergence and computational expenses 

The system of equations was solved numerically using modified 
version of FNAS2D code developed at CIAM [42] (description of the 
code is given in Sec.B.3 of Appendix B) . 

Details of calculations were basically the same as those 
reported in [27]. Only one additional modification was applied. 
The fully coupled solution procedure was used to increase 
convergence rate of the QL approach calculations compared to that 
reported in [27] (see detailed in Sec.B.3 of Appendix B). It 
allowed, by order, decrease required number of iterations and 
increase Courant number for QL calculations compare to the 
calculations reported in [27]. 

Main serie of calculations was done using nonuniform grid with 

90 cells in transversal direction and 100 cells in longitudinal 

direction (Fig. 27). Grid was clustered to the lower wall and to 

backward - facing step in accordance with geometrical 

progressions. The progression factors were chosen automatically to 

provide the given sum of progression but they were lower than 1 . 1 . 

In initial cross-section, minimal cell dimension near the lower 

wall was Ahmin =10” 5 m (Ahmin /h=2 . 5*10~ 2 ) . For the exit 

y y 

cross-section, minimal cell dimension in vertical direction was 
chosen to be equal to 2*10~ 5 m. This provided value y + =yU /v ~3 . 2 

J 2 Z w 
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for the nearest to the lower wall computational point (see 
Fig. 30). More detailed description of the grid generation 
procedure was given in [27] . 

The calculations were done with Courant number up to 400 for 
the FL approach and 1500 for the QL one. The convergence was 
estimated by the L 2 norm for the residual of continuity equation. 
The L 2 norm behavior vs iteration number is given in Fig. 31. In 
the case of FL approach calculations, this norm decreased 3.5 
orders during 600 iterations. In the case of QL approach 
calculations the L 2 norm dropped »3 orders after 200 iterations. 

The computations were performed using workstation HP 9000/735. 
The CPU time requirements are given in Table 5. 


Table 5. Computational expences for flowfield calculations 



QL approach 

(21 reactions, 
11 species) 

FL approach 

(21 reactions, 
11 species) 

Grid 

90x100 

90x100 

CPU time per grid cell per 
iteration 

0 . 007sec 

0 . 0013sec 

Total number of iterations 

200 

600 

Total CPU time of computation 

~3 . 5 hours 

~2 hours 


Accuracy of computations 

We were not able to control accuracy of computations by simple 
increasing of the computational cells number due to limitations in 
operational memory of our HP work station. So, to provide such 
analysis the following series of calculations was performed. 

The computational domain was divided by subregions and 
calculations with the patched grids were performed to estimate 
influence of the discretization in longitudinal direction. The 
final grid with 90x300 points was designed to provide very fine 
grid near the injector lip (100 points in the longitudinal 
direction) . Two another regions were joined with this first region 
in such a manner to provide the total number of grid points in 
longitudinal direction to be equal to 300. 

Additional calculations were performed also to estimate the 
sensitivity of the results to the grid accuracy in the lateral 
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direction. For this purpose, grid was adapted in such manner to 
increase the number of grid points in the mixing layer with 
something more rough grid in the near wall region. The grid 
adaptation to the stoichiometric line was realized. 

The obtained influence of the grid on the species mole 
fractions in cross-section x=0.356m is illustrated by Fig. 32a for 
FL approach, and, by Fig. 32b - for QL approach . It is seen that 
discretization influence on the results of calculations is «5%. 
There is no any qualitative difference between "main" (90x100) and 
"fine" (90x300) grid calculations. The mostly significant 
influence is confined in some shift of the QL predictions to the 
air-side of the mixing layer. It should be mentioned also that the 
main grid influence is connected with discretization in 
x-direction (especially for QL approach) . 

Of course, these estimations are slightly ambiguous. However, 
they provide to feel the influence of the grid accuracy on the 
results . 

RESULTS OF COMPUTATIONS 

The obtained H 2 0 mass fraction contours are given in Fig. 33 
for both FL and QL approaches. It is seen that the flamelet model 
predicted self - ignition point location 0.08m downstream hydrogen 
injector. The self - ignition point predicted by the QL approach is 
located 0.12m downstream the injector which is more close to its 
location in experiments (x»0.18m). The significant difference in 
H 2 0 mass fraction distributions predicted by the QL and FL 

approaches downstream the ignition point is seen also. It is 
necessary to note that the peak water concentration value 

predicted by FL approach is 17% lower than that predicted by QL 
approach . 

The comparison of the obtained reactive species mole fraction 

distributions with the experimental data is given in Fig. 34 for 

test section x=0.356m where the composition measurements were done 

by Burrows and Kurkov. Unsatisfactory correlation between FL 

approach predictions and experimental data for H 0 and 0 mole 

2 2 

fractions is seen. For example, the »35% underprediction of the 
H 2 0 peak value is observed. The results of the QL approach 
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predictions correlates with the experimental data reasonably well. 
The difference between measured and predicted by QL approach H 2 0 
concentration peak values is lower than 20%. 

The significant underprediction of the water origin in the 
flame by FL approach leads to approximately the same («35%) 
underprediction of the mean static temperature peak value as it is 
shown in Fig. 35 (solid line). It is seen that QL approach 
calculations gave satisfactory result again (denoted by dashed 
line in Fig. 35) . 

Unfortunately, one can conclude that the flamelet approach 
gave unsatisfactory results for the Burrows -Kurkov test case. 

Roots of the obtained discrepancy were analyzed in [27] . It 
was found that the mostly probable cause of discrepancy was too 
high level of the mixture fraction variance cr 2 predicted by 
eq.(II.lO), which was resulted in overprediction of the turbulent 
fluctuations intensity INT=<r/z level inside mixing layer for the 
considering wall-jet configuration. The main reasons leading to 
such conclusion are presented below. 

Let us start consideration from the discussion of the 
chemistry/pdf relative role in total FL approach predictions 
budget for the Burrows -Kurkov test case. It is given in Fig. 36. 
Here the dashed line denotes the equilibrium chemistry solution 
for water mole fraction * (eq) plotted vs mixture fraction z. It 
was obtained based on the assumption that both chemical 
nonequilibriumness and scalar field fluctuations are absent in the 
flow so it is the upper limit for possible water concentration 
distributions. The obtained in cross section x=0.356m 
instantaneous flamelet model solution for water mole fraction 
^H 2 o ' ( z < nS ' P s ) is plotted vs mixture fraction z by fine solid 
line. The averaged distribution of water mole fraction: 


I (fl ) , , , 

*H20 P (2)dZ 


is plotted in Fig. 36 by fat solid lines. 

It is seen approximately 50%-50% input of both chemistry 
nonequilibriumness and averaging procedure into the resulting 
averaged water mole fraction distribution X^^ predicted by the 
FL approach. Here averaged concentration peak value is 35% lower 
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than the peak value for distribution. Such a prediction 
contradicts to experimental data of Burrows and Kurkov denoted in 
Fig. 36 by crosses. It is seen that experimentally measured H 2 0 
mole fraction distribution is very close to the equilibrium limit 
*H 2 o’ This m eans that the role of both chemistry 
nonequilibriumness and scalar field fluctuations is very low for 
the considering test section of Burrow-Kurkov experiment. At the 
same time the flamelet approach overpredicted role of both 
effects. So one could expect that the flamelet model failure was 
due to failure of the used detailed kinetics model and/or due to 
failure of the averaging procedure. 

The additional tests demonstrated that version concerning 
detailed chemistry approximation failure for considering test 
conditions is improbable. The reasons here are as follows. The 
results of the flamelet model calculations varied very weakly when 
the Miller-Bowman chemistry approximation [32] was substituted by 
the Warnatz [34] one. At the same time the averaging of the 
equilibrium chemistry distribution * (eq) (using obtained in 
calculations distributions of z and a 2 ) gave the mean 


*) 

This was done as follows. The measured in Burrows -Kurkov 
experiment distributions of N 2 mole fraction was used to re-plot 

reactive species mean mole fractions vs mean mixture fraction z 
instead of space coordinate y. It was expected that the molecular 
nitrogen N 2 does not react with other species. Hence its mass 

fraction C n2 has to be approximately linear depending on mixture 

fraction: C N2 =C* 2 (l-z) if difference in molecular dif fusivities is 

neglected. Here C* 2 is the molecular nitrogen mass fraction in the 

air flow. The fluctuations can not change this dependence due to 
its linearity and it remains valid for mean values C =C A (1-z) 

N2 N2 

We used this relation to obtain the mean mixture fraction 

distributions z (y) in the test sections. We re-plotted 

experimental data of Burrows -Kurkov in a form of x a on z, using 

the experimentally measured distributions of the reactive species 
mole fractions * a (y) and the obtained relation z (y) . Such a 

re-plotted H 2 0 mole fraction distribution (denoted by crosses) is 

given in Fig. 36. The experimental N z mole fraction distribution 

had non-monotonic part (three measured points) at the fuel -lean 
edge of the mixing layer (please see Fig. 34) which can be referred 
to instrumental errors. That is why these 3 points were rejected 
from the analysis. 
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distribution which was very close to the averaged flamelet 

model solution * <FL) . 

H20 

One could expect improper thermochemistry model operation due 
to using of the simplified correlation U=U(z) which could 
introduce some errors in enthalpy distribution. However we 
controlled results sensitivity to U=U(z) correlation influence. 
This influence was very weak (please see Figs.26a-d) . 

So one can expect that the roots of the failure are 
connected with the averaging procedure. 

First of all we have verified role of the shape of our pdf 
model . For this purpose we have performed averaging of the 
instantaneous solution using two another approximations for p(z). 
The first one was widely used ^-distribution [43] . As the second 
approximation for pdf we have used analytical solution of the 
mixture fraction pdf equation obtained in [44] for variable 
density homogeneous turbulent flows. The obtained averaged 
distributions were slightly different but minor difference for the 
peak values of averaged water concentration was found (Fig. 37) . 

The "smoothing" of the species and temperature instantaneous 
distributions by probability density function is governed by the 
parameter INT=<x/z. The results of averaging are mostly sensitive 
to INT level in near- stoichiometric regions (z~z ) where 

S 

significantly non-linear variation of the instantaneous reactive 
species distributions and temperature vs mixture fraction is 
observed. At the same time the main discrepancy between FL 
approach averaged distributions and experimental data was observed 
in near-stoichiometric regions also (see Fig. 36) . This feature 
sent us to analyze obtained in calculations INT distribution. This 
distribution vs mean mixture fraction z is given in Fig. 38 for 
test section x=0.356m. It is seen that the eqs . (II . 9) , (II . 10) 
generated INT distribution with a level of about 150-200% in the 
near-stoichiometric regions (z^zj . This prediction of 
"concentration subsystem" (II. 9) , (II. 10) of the turbulence model 
is quite questionable for the considering planar wall jet 
configuration. Results of presented analysis of Burrows -Kurkov 
experimental data allows to expect that INT level in their 
experiments was sufficiently lower. The reasons here are as 
follows. It is known that approximately the same INT level 
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(*200%) , as it was predicted by eqs . (II . 9) , (II . 10) , is typical for 
unconfined flames in still air where it provides significant 
difference between instantaneous and averaged distributions of 
parameters [1] . For example, it provides *400K (*20-25%) 
difference between instantaneous and averaged peak temperature 
values in laboratory H 2 /air round jet diffusion flames which was 
approved both by experiments and calculations [45] . 

The possible overprediction of INT level explains also why the 
instantaneous flamelet solution differed so significantly from 
equilibrium limit (see Fig. 36). The nonequilibriumness for FL 
approach is governed by approximation (1.15) for conditionally 
averaged value of scalar dissipation at the stoichiometric 
surface : 

N s = N | z = zs 

t , ~ 

7 | Z=2s 

_ s g 

The value of N t ~ INT where exponent B is varied in the range 2-4 
i.e. N^ value is very sensitive to the overprediction of 
fluctuations intensity level. So the overprediction of the 
fluctuations intensity leads to overprediction of the N s value and 
as a result to overprediction of chemical nonequilibriumness role 
by flamelet model equations (1.5) . 

We estimated possible range of the intensity reduction which 
can improve the results of the FL approach predictions. The 
INT =INT(z) distribution obtained in FL approach calculations was 
diminished "by hand" by a factor i(t = 0.5 and 0.25 sequentially 
(INT=^- INT° ) . The influence of INT decreasing on the predictions 
of water mole fraction values is illustrated by Fig. 39. It is seen 
that improvement can be obtained if intensity of fluctuations 
would be *3-4 times lower than that obtained in current FL 
approach calculations. 

Based on the results of the presented here analysis, we expect 
that the fluctuations intensity in the conditions of 
Burrows -Kurkov experiment was low enough (especially in 
near- stoichiometric and fuel -lean regions) and subsystem 
(II . 9) , (II . 10) failed to predict it accurately. We expect also 
that the low level of mixture fraction fluctuations provided 
relative success of the QL approach obtained not only in our but 
in many other QL approach calculations for the conditions of 
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Burrows -Kurkov test case (for example [4,5,46]). 

We expect that the overprediction of the I NT level in our 
calculations was mostly due to overprediction of mixture fraction 
variance <r 2 level by semi-empirical eq.(II.lO). Previously, this 
equation was verified mostly for round jet configurations and for 
low Mach numbers. We expect that performed tests demonstrated that 
this equation requires additional verification and improvement 
for complex flow configurations (compressible flows, wall jets, 
etc.). We expect also that adopted approach for the calculation 
of the intermittency factor ? requires additional re-examination 
and (may be) improvement. The intermittency factor j was 
calculated from the simplified relation 2r=min{l;A/ (1 + INT 2 ) } (see 
eq.(1.12)) which was previously obtained and verified in [1,29,30] 
for incompressible flows only. At the same time some of the 
available experimental studies performed for compressible flows 
demonstrated decreasing of the intermittency effects role with the 
increasing of the Mach number due to the role of compressibility 
[47,48] . We hope that such improvement of the "passive scalar 
statistics" block in our turbulence model will allow to improve 
accuracy and reliability of FL approach predictions for reactive 
scalars and temperature without modification of flamelet 
equations (1.5) . 

However we have doubts concerning possibility to improve 

predictions of self -ignition point location by the existing FL 

approach version. This conclusion is based on the fact that 
existing FL approach considers only one possible cause responsible 
for ignition delay i.e. too high rate of reagents mixing in the 
vicinity of the injection point which limited rate of chemical 

reactions (due to reaction zone cooling) . At the same time there 
is another possible cause responsible for ignition time delay i.e. 
pure chemical kinetics limitations. Both experimental data and 
results of QL approach calculations demonstrate that such kind of 
self - ignition was observed in Burrows -Kurkov experiment. So the 
improvement of the flamelet approach for proper prediction of 

self - ignition point location in such situations requires further 
upgrading of the flamelet equations (see Sec. II. 3). 
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II. 3 COMPARATIVE ANALYSIS AND CONCLUSIONS 

Current Section summarizes our conclusions concerning 
performed CFD tests of the model. Computational efficiency of the 
FL approach, accuracy of its predictions, roots of discrepancy and 
our suggestions for the model improvement are discussed. We 
presented also estimations of the flamelet model capabilities in 
comparison with the another turbulent combustion models (Evolved 
PDF and Assumed PDF approaches) . The necessary data for 
comparative estimations were taken from study [12] where Evolved 
and Assumed PDF approaches were applied for calculations of the 
same test case [24] . 

COMPUTATIONAL EXPENSES 

The summary of computational expenses required for FL approach 
calculations is presented in Fig. 40. Data were taken from Tables 
2-5 of Sec . II . 1, II . 2 . Additionally it was taken into account that 
PC AT 4 8 6DX2 / 6 6MHz productivity is *1/10 of HP 9000/735 work 
station . 

It is seen that application of the flamelet model in CFD does 

not require any significant computational resources. For example, 

we found that generation of the sufficiently accurate (accuracy is 

about 0.3-0. 5%) flamelet library for the detailed H oxidation 

2 

chemistry (11 species, 21 reactions) required -20% of the total 
CPU time necessary for the flowfield calculations if PNS marching 
solver is used. The relative CPU time expenses for FL library 
productions decreased to the value -4% if full NS solver is used 
for flowfield calculations. 

The calculations of the test cases [24,25] demonstrated also 
that the FL approach usage allows to decrease the computational 
time per cell in -1.5 times for the PNS solver and in *5 times for 
the NS solver in comparison with QL approach. 

The flamelet libraries does not require any significant memory 
for storing (~l-3Mb) . 

The rough comparative estimation of CPU expenses for FL, 
Evolved PDF and Assumed PDF modeling is given in Fig. 41. Here the 
expenses for Evolved PDF are adopted as 100% and y-axis is 
logarithmically scaled. The Cray-YMP CPU requirements for Evolved 
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and Assumed PDF calculations of Beach test case were reported in 
study [12] . The CPU expenses for FL approach (600sec of HP 
9000/735 CPU; 50x200 grid) were obtained in FL calculations based 
on full NS solver FNAS2D (described in Sec.B.2 of Appendix B) of 
the Beach test case also. The minimal realistic assumption that 
productivity of Cray-YMP is only 10 times higher than productivity 
of our HP work station was adopted in comparative estimations. It 
is seen that computational expenses for FL approach by orders 
lover both Evolved and Assumed PDF modeling. 

Based on the reported results, one can conclude that performed 
tests demonstrated high computational efficiency of the FL 
approach . 

ACCURACY^ OF PREDICTIONS and SUGGESTIONS FOR MODEL IMPROVEMENT 

Unfortunately, results of the performed tests does not allow 
to make so definite conclusion concerning accuracy of the model 
predictions as it was made about its computational efficiency. 
Moreover, these tests demonstrated that model requires further 
improvement . 

The mostly important results of the tests are summarized in 
Fig. 42. Here results of our FL approach water concentration 
predictions (taken from Figs. 21 and 34 ) are plotted for both 
considered test cases. The QL approach (current), Evolved PDF 
modeling [10,12] and Assumed PDF modeling [12] predictions are 
plotted also. 

It is seen that correlation between the FL approach and the 
experimental data was quite different for the considered test 
cases . 

The FL approach gave satisfactory results for the Beach test 
case (Fig. 42a) . As a rule, predictions of species distributions 
were possible with accuracy better than *20%. Here FL approach 
predictions looks quite well compare to those obtained for Beach 
test case by another approaches . 

However FL approach predictions were unsatisfactory for the 
Burrows -Kurkov test case (Fig. 42b) . The serious discrepancy 
between experimental data and results of calculations was found 
i.e. the »35% underprediction of the H z O and temperature values; 
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large "smoothing" of 0 2 distribution, etc. It is seen that Evolved 
PDF modeling for the considering test case provided best result. 
The QL approach predictions were satisfactory. 

The analysis of the discrepancy roots was done in [27] and it 
was repeated in Sec. II. 2. Based on its results, we can conclude 
that the cause of discrepancy was not directly connected with the 
flamelet equations. We expect that it was due to overprediction of 
the mixture fraction intensity level by our turbulence model. The 
main doubts are addressed here to the accuracy of semi-empirical 
eq.(II.lO) for mixture fraction variance. We hope that even 
current version of FL approach can get satisfactory results for 
the predictions of the reactive species and temperature 
distributions after improvement of the "passive scalar statistics" 
block in our turbulence model. To obtain this goal, we intend to 
re-examine accuracy of closure approximations and role of 
neglected terms in (II. 10) using available experimental data on 
turbulent mixing for flows with and without burning. Based on this 
analysis, we expect to propose modifications of the model 
equations for mixture fraction variance and intermittency factor 
which will allow to increase reliability and accuracy of 

predictions for complex flow configurations like that which was 
considered in current study (compressible wall-jet duct flow) . 

The performed tests demonstrated also that existing version of 
the flamelet model does not allow to predict accurately 
combustible flows with large ignition delay distance. This is 

mostly due to improper model operation in the vicinity of 

self-ignition point since convective terms in the flamelet model 
equations were dropped out. This approach is approximately valid 
in the mixing region upstream ignition zone and in the 
"well -developed" combustion region downstream this zone. But it is 
violated in the "self -ignition" region, where reaction zone is 

thick and basic assumptions of the approach can violate. We expect 
that this disadvantage can be overcome by introducing of 
convective terms into the flamelet equations based on the same 
kind of reasoning as it was done in Condition Moment Closure (CMC) 
approach proposed by R.Bilger (Sydney University, Australia) [49] 
and A. Klimenko (CIAM, Russia) [50] . We expect to investigate 
capabilities of such model upgrading for predictions of flames 
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with significant ignition length delay. The role of the scalar 
dissipation fluctuations is additional important effect which can 
significantly influence flamelet model predictions in the vicinity 
of the ignition point. The role of this effect was neglected in 
current study. So we expect to introduce pdf model for scalar 
dissipation into our averaging procedure and to investigate 
influence of the scalar dissipation fluctuations on the 
predictions of the ignition point location and averaged 
distributions of reactive species in vicinity of this point. 

The discussed here two directions (improvement of the passive 
scalar statistics modeling and upgrading of the FL approach for 
flames with large ignition delay) were proposed to National 
Aeronautics and Space Administration as main tasks for the next 
step research directed toward further improvement of the FL 
approach . 
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Part III. SUPPLEMENTED STUDIES 


III.l SIMPLIFICATION OF THE FLAMELET CALCULATIONS BASED ON PARTIAL 
EQUILIBRIUM CHEMISTRY ASSUMPTION 

Detailed kinetics of hydrogen oxidation consists of two main 
subsystems. The first one couples relatively fast two-body chain 
branching and chain-propagation reactions of radicals production: 


H + OH = H 0 + H 

2 2 

OH + OH = H 0 + 0 

2 


tvo-body reactions 


(III.l) 


The second subsystem couples three-body radical recombination 
reactions which determine reaching of the chemical equilibrium by 
the reacting system : 


H + OH + M = H 0 + M 

2 

H + 0 + M = OH + M 


three-body reactions 


(III .2) 


Here third body M is required to carry out away liberated energy. 

In many classes of high-enthalpy reacting flows, induction 
times for radicals production via two-body reactions are short 
compare to their consumption in three-body recombination reactions 
[59,60]. This feature of combustion chemistry was used previously 
to develop Partial Equilibrium approach (PEq) for simplification 
of the detailed kinetics approximation [61,62] . 

PEq is based on the assumption that two-body subsystem (III.l) 
is equilibrated while all departure from the chemical equilibrium 
is governed by three -body subsystem (III. 2). The equilibrium 
conditions for two-body reactions allows to reduce number of 
required scalars for description of combustion chemistry to only 
two. Usually, mixture fraction z and primary kinetic variable C* 
are used as such variables [63] . 

Range of PEq applicability is limited by the moderate 
deviations from the equilibrium. Significant increasing of the 
chemical nonequilibriumness leads to PEq violation due to 
increasing of the induction times of two-body reactions. However, 
benefits of PEq two-variables formalism were used successfully for 
the simplification of chemistry calculations in conventional 
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turbulence modeling of both sub- [43,63] and supersonic [64] jet 
flames . 

In current study, we tried to apply PEq two-variables 
formalism for simplification of the flamelet calculations. It was 
possible to propose simplified semi-analytical method for 

approximate solution of the flamelet equations at moderate 
deviations from the chemical equilibrium. Description of the 
simplified method and account about its accuracy are presented 
below . 

Simplified procedure for FL equations solution 

Followed by [61,63], it is assumed that thermochemical system 
for H 2 oxidation can be approximated using six reactive species 
(H 2 , 0 2 , H 2 0, 0, H, OH) consisted of two elements (H,0). The 
molecular nitrogen N 2 is considered as passive contaminant. 

The PEq assumption about equilibrium of two-body reactions is 
used. Due to 0 and H atoms conservation laws, it is possible to 
select only three independent stoichiometric equations between six 
reactive scalars from the two-body reactions subsystem. The 
following set was selected: 

H + 0 = OH + 0 

2 

H 2 + OH = H 2 0 + H (III. 3) 

OH + OH = H 0 + 0 

2 

It is expected that all deviations from chemical equilibrium 
are governed by the three-body recombination reactions: 

H + OH + M = H 0 + M 

2 

H + 0 + M = OH + M ( 1 1 1 . 4 ) 

H + H + M = H + M 

2 

For further consideration, let us denote H , 0 , HO, 0, H, OH 

2 2 2 

mass fractions as C , . . .,C respectively. 

Let us introduce primary kinetic variable C* as linear 
combination of reactive species so as it to be changed only by the 
slow recombination reactions (III. 4). The choice of C* is based on 
method proposed in [62]. It is as follows. 

Flamelet equation for a-th specie from system (1.5) : 

2 

N s d_|^ + r (p,T f C , ...,C ) =0 a=l , . . . , J=6 (III. 5) 

dz 
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is multiplied by arbitrary value A ft and the resulting relations 
are added. As a result one has: 

i 2 / J \ J 




dz‘ 


f Ya .cl = - Y A -R 
I L a a) L a a 

v a = i ’ a = i 


(III. 6) 


Using the relation: 


n 

R = U Y (i> b - v f ) -RR 
a a L am am m 

m = 1 

it is possible to transform right hand side of eq. (III. 6) as 
follows : 


J 

M J 





y a -r = 

L a a 

l rr . l 

II b 

( (v 

OCm 

- V { ) (1 A ) 
am *a a 


(III .7) 

oc = \ 

m = 1 a= 

1 



where RR is the 

m 

total 

rate (difference 

between 

forward and 

backward rates) of 

m-th 

elementary reaction 

from set 

of (III. 3), 


(III. 4); ^m'cL" stoichiometric coefficients for a-th species in 
m-th elementary reaction; u a is molecular weight of a-th specie. 

If the coefficients at the bimolecular reaction rates in 
(III. 7) to equate with zero, one can obtain three linearly 
independent equations: 

j 

0 (III. 8) 


Y ( (l/ - v* ) U A 
L am am a a 


a = i 


with regard to value of . It is convenient to assume that 
coefficient A i at ^ (H 2 concentration) is equal to unity and 
coefficients A 2 ,A 3 at C 2 ,C 3 (C> 2 and H 2 0 concentrations) are equal 
to zero. Then from solution of the equations (III. 8) it follows: 


C = C + 
i 




h. 


C 4 + T 






C S + T 


M. 


h. 


(III. 9) 


4 *""5 " '“6 

By the definition (III. 9), primary kinetic variable C* obeys 
equation : 

* 


N 


d 2 C* 

dz 2 


+ R (p,T,C lf . . . ,Cj) = 0 


Here, effective oxidation rate R* has the form: 


* 

R = 


2/l l (RR + RR + RR ) 

1 1 2 3 


(III . 10] 


(III. 11) 


where RR x , RR 2 < RR 3 _ total rates of three-molecular reactions 
(III. 4); subscripts 1-3 in the values of RR^ correspond to the 
ordinal number of the reaction in (III. 4). 
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Generally speaking value of the effective oxidation rate R* is 
expressed through p and T using Arrenius approximation 
of rates for reactions (III. 4). So, the following set of equations 
should be used together with eq.(III.lO) to obtain solution: 


i ) atoms conservation lavs in a form: 
(H atoms) =* 


(0 atoms) 


(ZH f -ZH a ) -z+ZH a = C + l-C + C + -i--C 

1 9 3 5 17 6 

(Z0 F -Z0 A ) -Z + Z0 A = C + s'C + C + -c 

2 9 3 4 17 6 


(III . 12) 


ii) flamelet form (1.5) of energy conservation equation: 


H=(H F -H A )z + H A 


(III . 13) 

Hi) equilibrium conditions for two-body reactions (III. 3): 


n (^H 

a=i v a } 


{v 


Otm 


V* ) 
0 dm 


eq 


= K 


iv ) thermal equation of state: 


P= PET l C a / Ua 


(hi . 14 : 


(III .15) 


a=i 


v) definition of the primary kinetic variable C 


C*= C + 

l 

Here : ZH F , 
fuel flow 


3 

— — C + — 
u 4 2 




C S + T 




(Ill . 16) 


Z0 are hydrogen and oxygen atoms mass fractions in the 
(z = l) ; ZH? Z0 A are hydrogen and oxygen atoms mass 


eq 

fractions in the oxidizer flow (z=0) ; H is total enthalpy; K 

m 

equilibrium constant; R is universal gas constant. Additionally, 
similarity between mixture fraction conservation equation and 
atoms conservation equations is taken into account in (III. 12). 

Necessity to solve large subsystem of additional algebraic 
equations together with differential eq.(III.lO) decreases 
computational benefits of the PEq approach. So, based on results 
of [20] , the following simplified procedure was used. 

The effective oxidation rate R* is approximated by a simple 
formula : 


R = - k- (C*-Ce) n , (III. 17) 

where k and n are some constants, C* is equilibrium value of C*. 

The factor k and exponent n for approximation (III. 17) can be 
found from the algebraic subsystem ( III . 12) - ( III . 16) . This 
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subsystem is solved parametrically for various (and given) values 
* # * 

of AC =C -Ce. As the result, reactive species, temperature and 

density are found as the function of mixture fraction z, increment 
... * 

of the primary kinetic variable AC , and pressure in the reaction 
zone Ps : 


C =C (z, AC , P 
a a s 


T=T ( z , AC , P 


(III . 18) 


The obtained relations are substituted into the relation (III. 11) 
for the effective oxidation rate R*. As the result, one has: 

R*=R* ( z , AC* , P ) (III. 19) 

s 


Further k and n are obtained by approximation of (III. 19) . 

Using approximation (III. 17), equation (III. 10) for primary 
kinetic variable is written as: 
d 2 C 

N s — - k (C - Ce ) n = 0 (III. 20) 

dz 

It is seen that this equation can be solved independently on 
algebraic subsystem ( III . 12 ) - ( III . 16) in parametric form on N s . 

Two approaches were formulated previously for approximate 
analytical solution of eq.(III.20). The first one was reported in 
study [20] . It is based on the expansion over small parameter 
e=N s /k<<l which is contained in eq.(III.20). Here, for the first 
approximation, solution has the form: 



if d 2 C * 
k dz 2 



(III. 21) 


The second analytical approach for solution of (III. 20) was 
proposed previously in our group [1] . Here, the equilibrium value 
of primary kinetic variable Ce is approximated by formula: 


Ce = k(z-z)0(z-z) 

Os s 

where k =(l-z ) 1 and 6 is Heaviside step function. This 

US 

approximation allows to obtain the following analytical solution 
in the vicinity of the stoichiometric surface: 


C=k (z-z )0( z-z ) + fi- 

Os s 


abs (z-zs) 


where Damkohler number is defined as Da= N s k ( 1 n) /~k.; Q and * are 
some constants depended on n. 
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Unfortunately, tests, performed in the course of the current 
study, demonstrated that both solutions have very limiting range 
of applicability restricted by extra-low departure from chemical 
equilibrium (N s <0 . Olsec -1 ) . So, they are useless for practical 
calculations . 

Keeping in mind this finding, we used approach where boundary 
problem for equation (III. 20) was solved numerically. Further, 
substituting solutions for C*, we re-calculated reactive species 
concentrations and temperature vs z and N s using solutions 
(III. 18) of the algebraic subsystem (III . 12) - (III . 16) . 

So, simplified PEq method consisted of three particular steps: 

i) At the first step, algebraic subsystem (III . 12) - (III . 16) was 

solved parametrically for different values of primary kinetic 

variable increment AC*=C*-C*. The concentrations and temperature 

were obtained in a form (III. 18): 

C =C (z, AC* , P ) 
a a ' ' s 

T=T ( z , AC* , P ) 

S 

These dependencies were used to obtain factor k and exponent n for 
approximation (III. 17). 

ii) At the second step, single differential equation (III. 20) was 

solved numerically for various values of scalar dissipation N s . 

Scalar sweeping was used for solution of the finite difference 

equation. As the result, increment of primary kinetic variable 
* * * 

AC =C -Ce was obtained in parametric form on mixture fraction z 
and N s : 

AC*= AC* { z , N s , P s , BC) (III. 23) 

where BC denotes boundary conditions at z=0 and z=l ends. 

Hi) At the third step, results of parametric calculations 
(III. 18) for algebraic subsystem which were obtained at the first 
step and solution (III. 23) were used together to re-calculate 
dependencies for reactive species and temperature vs z and N s in 
conventional for flamelet approach form (1.9): 

C =C ( z , N s , P , BC) 

OC OC s 

T=T ( z , N s , P , BC) 

S 

So, single ordinary differential equation was needed to solve 
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in simplified PEq method instead of full system of flamelet model 
equations (1.5) and without internal stiffness in the source term. 

Examples of calculations and accuracy 

We obtained flamelet solutions using simplified PEq method for 
conditions of Beach and Burrows -Kurkov test cases [24,25] 
considered in Part II of the Report. The kinetics constants for 
recombination reactions (III. 4) were taken as those in original 
Miller-Bowman detailed kinetics scheme presented in Appendix C. 
The equilibrium constants for two-body reactions (III. 3) were 
taken in accordance with [33] . 

The calculated by PEq method best fit approximations for 
effective oxidation rate R were close one to another for both 
considered test cases. It provided approximation for R* as: 

* s * * ? 

R = - 5.5-10 - (C -Ce) 

Using obtained approximation for R*, we calculated reactive 
species and temperature distributions vs z for different values of 
scalar dissipation N s in the range 0<N s =a00sec _1 . 

The accuracy and range of applicability for PEq method were 
controlled by comparison with "accurate" flamelet solutions 
obtained based on detailed kinetics model. The accurate solutions 
of the flamelet equations were obtained using Miller-Bowman 
detailed chemistry model [32] (see Appendix C) without assumptions 
concerning equilibria of the two-body reactions*'. These solutions 
were used to calculate effective oxidation rate R* and primary 
kinetic variable C directly from the definitions (III. 9), 
(III. 11) and without any assumptions adopted in PEq method. 

The accuracy of the simplified PEq method for calculation of 

the effective oxidation rate R* is illustrated by Fig. 43. Here, 

both approximation for R* provided by simplified PEq method (blue 

line) and results of accurate flamelet calculations based on 

detailed kinetics and definition (III. 11) (black points in the 

rose area) are presented for the whole range of z and N s 

variation. It is seen that approximation (III. 17) allows to 
* 

calculate R with the accuracy about 10% in the range of N s 

*) 

In reality, these accurate solutions were taken from the 
flamelet libraries generated for FL approach CFD tests reported in 
Part II. 
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variation 0<N s ^100sec 1 . Further increasing of N s leads to rapid 
loosing of accuracy by approximation (III. 17). 

The primary kinetic variable increment AC* can be calculated 
using PEq method with the accuracy within 6% in comparison with 
accurate flamelet solutions in the range 0<N s il00sec _1 as it is 
shown in Fig. 44. The aforementioned range of N s variation 
corresponds to AC variation in the range 0<AC*^0.02 

The comparison of the simplified PEq method predictions for 
species mass fractions and temperature with those provided by 
accurate solution of the flamelet equations based on full detailed 
kinetics scheme are presented in Figs.45a-d (Beach test case) and 
46a-c ( Burrows -Kurkov case) . The results of calculations provided 

by simplified PEq method are plotted by dashed lines. The results 
of accurate flamelet calculations are presented by solid lines. It 
is seen that simplified PEq method allows to calculate with good 
accuracy reactive species mass fractions and temperature at 
moderate deviation from equilibrium. The accuracy of PEq method 
was better than 1.5% for stable substances, better that 3% for 
temperature, and better than 10% for main radicals. The obtained 
inaccuracy of PEq method is small enough. It is close to 

scattering of results obtained based on different detailed 
kinetics schemes as it was reported in Part II. It was possible 
to obtained flamelet solutions by simplified PEq method with such 
accuracy in the range of scalar dissipation variation 

0<N s ^100sec _1 . The accuracy provided by the PEq method decreased 
for higher values of N s . Of course it can not be used for 

calculations at highest values of N s in extinction/ignition region 
where PEq approximation became worth. However, performed 
calculations shows that more than 50-60% of the FL libraries can 
be generated based on outlined simplified PEq method. Its range of 

applicability by four orders higher than that provided by 

analytical solutions ( III . 21) , ( III . 22) . It is interesting to note 
that approximation (III. 17) for effective oxidation rate has 
sufficiently wide range of applicability (in current calculations 
up to N s /N *0.15, where N is critical value of scalar 

cr cr 

dissipation) . So, may be, such kind of analytical approximations 
could be useful for simplification of the terms contained 
derivatives over reactive scalars in Evolved PDF modeling. 
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Ill .2 FLAMELET PREDICTIONS FOR EXTINCTION/ IGNITION PHENOMENA 

As a rule, combustion chambers and burners operate at moderate 
values of the incoming flow temperature (Ta«450-800K) . It is 
well-known [65] , that it is impossible to provide stable 
combustion regimes without implementation of flameholding systems 
at such conditions. Induction times for reaction chemistry are too 
high and reactions acceleration requires increasing of incoming 
flow temperature by special devices. Flameholders (bluff bodies, 
flow swirling, etc.) are used to transfer some heat from zones 
with "well -developed" combustion back to the fresh mixture. 
Loosely speaking such a transfer helps to change the "state of 
mixing" to the "state of burning". 

The increasing of the incoming flow temperature decreases 
rapidly induction times of chemical reactions. There are some 
boundary conditions of operation where flameholding became 
useless. The combustion wave can self-develop at such conditions 
without additional backward heat transfer and there is no need to 
use flameholders. For example, such regimes could be expected in 
scramjet combustors due to high value of the incoming static 
temperatures . 

We tried to estimate location of the boundary between 
aforementioned combustion regimes using flamelet model. For this 
purpose we considered features of the flamelet model 
ignition/extinction predictions for H /air flames in different 
ranges of operational conditions . Flamelet model calculations were 
performed based on Miller-Bowman detailed kinetics model for H 

2 

oxidation and thermodynamics approximations for species 
enthalpies presented in Appendix C. 

Typical dependencies of the maximum temperature value in 
H 2 /air diffusion flame T m vs scalar dissipation N s , which were 
obtained from the flamelet model calculations, are presented in 
Fig. 47 together with the regime conditions. 

It is seen (CASE A in Fig. 47), that flamelet eqs.(1.5) have 
two solutions at moderate values of the incoming temperatures of 
reactants. The first one corresponds to the burning regime (red 
line) . The second solution corresponds to the mixing of reactants 
accompanied by the very slow oxidation (blue line) . It is seen 
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also that there are three different regions which are separated 
one from another by parameters N (1) and N (2> . We will refer to 

cr cr 

these parameters as critical values of scalar dissipation at flame 
extinction (N^ 1 ) and flame ignition (N^ 2, )* ) . Only mixing 
solution exists at very high values of scalar dissipation in the 
reaction zone N s (region N s >*r* ) due to too high rate of reagents 
mixing compare to their consumption in chemical reactions. Both 
burning and mixing solutions can be found if value of N s obeys 
condition N <2) <N s <N ll) . Only burning solution exists if N<N (2) . 

cr cr 

So, for this combustion regimes, transition from the mixing to the 
burning requires flameholders "to jump" from mixing solution to 
the burning one and to overcome finite chemistry limitations in 
fresh mixture. We denoted such combustion mode as CASE A. 

Qualitatively another results were found at high values of the 
incoming reagents temperatures (CASE B in Fig. 47). It is seen that 
flamelet solution is unique for such conditions at any values of 
the scalar dissipation in the reaction zone N s . So, transition 
from mixing to the burning regime occur continuously and it does 
not require flameholders since chemistry is rapid enough to 
provide self -ignition . We denoted such combustion mode as CASE B. 

It is seen that boundary which separates CASE A from CASE B 
combustion regimes can be introduced by the condition: 


n <2) =n (1) 

cr cr 


(III. 24) 


Both critical values of the scalar dissipation depend on incoming 
fuel and air temperatures (T ,T), pressure in reaction zone P 

fa s 

and detailed chemistry of fuel oxidation (kind of fuel) . These 
dependencies can be found from the flamelet model calculations as : 


N (1) =N (1> (P , T , T , kind of fuel) 

cr cr s f a 

( 2 ) ( 2 ) 

U =N (P , T , T , kind of fuel) 

cr cr s f a 


(111.25! 


Substituting solutions (III. 25) into condition (III. 24) one has 
relation between thermodynamics parameters at the boundary which 
separates CASE A combustion mode from CASE B one in a form: 


T ,b, = *,(T <b :p' 6 -) 

a £ s 


(III . 26! 


* J We used single value of Ncr in CFD tests of the flamelet model 
reported in Part II. It corresponded to Nir value. 
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We estimated location of this boundary for the H 2 /air 
diffusion flames. Parametric flamelet calculations were done in 
the range of conditions summarized in Table 6 to cover possible 
range of parameters variation in scramjet combustors at different 
flight conditions. The air composition was taken as: C> 2 mass 
fraction 0.231 and N 2 mass fraction 0.769. 

Table 6. Initial data for FL calculations 


Ps ,MPa 

Ta , K 

Tf , K 

0.05-0.2 

300-2600 

300-1000 


We calculated dependencies (III. 25) for both critical values 
of scalar dissipation. The regimes, where values of N <2) and N (1) 

cr cr 

differed one from another less than 1%, were adopted to adjust 
location of the boundary between CASE A and CASE B combustion 
modes . 

Examples of the obtained parametric dependencies of N (1) and 
(2) cr 
N cr are given in Fig. 48a, b. It is seen (Fig. 48a) that the mostly 

important parameters governed transition from CASE A to CASE B 

regimes are the air and fuel incoming temperatures T , T . Their 

increasing rapidly decreases difference between N (1) (red lines) 
(2) . cr 
and N (blue lines) . More weak influence of the pressure P on 

the transition boundary was found in calculations (Fig. 48b) for 

the considered range of P variation (P =0 . 05-0 . 2MPa) . 

S S 

Summary of the obtained results is presented in Fig. 49. Here 
the obtained incoming air boundary temperature Ti b) which 
corresponds to the transition between combustion modes is plotted 
vs fuel temperature T^ b) for different pressures in the reaction 
zone . 

We compared obtained results with some typical scramjet 

combustor operational conditions at different flight Mach numbers 

to estimate, roughly, their possible practical significance. The 

summary of used scramjet operational conditions is presented in 

Table 7. They correspond to "hypothetical" constant dynamic 

pressure scramjet trajectory q =const~50-75KPa, where q is 

00 ^00 

dynamic pressure of free air stream. 
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Table 7. Expected scramjet operational conditions 
(flight trajectory q ro «50-75KPa) 


M 

00 

8 

10 

12 

15 

M 

air 

-3-3.5 

-4 

-4.5 

-5-5.5 

T , K 

air 

900-1000 

1100-1200 

-1500 

-1900 

P , MPa 

air 

0.08-0.12 

0 .1-0 . 07 

0.08-0.06 

0.06-0.05 


-3-3.5 


-1000 


-0.09 

-0 . 11 

-0.15 

-0 . 18 


Here: M is flight Mach number; M , T , P are the Mach 

00 air air air 

number, static temperature and static pressure for the air flow at 
the entrance of combustor; M , T* , P are the Mach number, 

H2 H2 H2 

total temperature and static pressure for the entering hydrogen 
flow . 

Results of comparative estimations are presented in Fig. 50. It 

is seen that they does not give hope. Flamelet model predicts that 

the CASE A combustion mode dominates up to the very high flight 

Mach numbers (M^-14) . Transition from CASE A to CASE B combustion 

regimes can be expected at near-orbital vehicle velocities 

(M -15-16 and higher) . 

00 

Of course, presented estimations are rough enough. They does 
not take into account accurately many important flow features 
(shocks, induction delay of the chemistry, etc.). Thus, this 
finding demonstrates that, now, performed calculations have mostly 
academic interest for the study of the flamelet solutions features 
for flame critical regimes rather than any practical application. 
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Ill .3 CONCLUSIONS ABOUT RESULTS OF SUPPLEMENTED STUDIES 

Two supplemented studies, related with the investigation of 
the flamelet model capabilities were done in the course of current 
research work . They were : 

° Consideration of the possible simplifications in flamelet 
calculations based on partial equilibrium assumption for the 
detailed oxidation chemistry . 

° Study of the FL predictions features for ignition/ extinct ion 
phenomena in high-enthalpy flovs. 

Performed investigation demonstrated that: 

1. It was possible to formulate simplified procedure for solution 
of the flamelet model equations at moderate (N s /N <0 . 1-0. 15) 

cr 

deviations from chemical equilibrium using Partial Equilibrium 
(PEq) assumption for the H 2 combustion chemistry. The simplified 
PEq method is based on numerical solution of single flamelet 
equation for primary kinetic variable C* where effective oxidation 
rate is analytically approximated using equilibrium conditions for 
two-body reactions, energy and atoms conservations laws [20] . 
Tests of its accuracy were done. It was found that accuracy of PEq 
method was better than 1.5% for stable substances, better that 3% 
for temperature, and better than 10% for main radicals in the 
aforementioned range of N s /Ncr variation. This finding allows to 
apply simplified PEq method for calculations of about 50-60% of 
the total FL libraries. We expect that this results are 
encouraging and it would be useful to test capabilities of PEq 
method in more wide range of regime parameters and for another 
fuels than that considered in current study. 

2. Feature of flamelet solutions in the vicinity of the flame 
ignition/extinction regimes of H 2 /air diffusion flames were 
studied . It was found two different combustion modes depending on 
the range of regime parameters. The first combustion mode (CASE A 
mode) corresponded to the regimes where flameholding is needed for 
stabilization of the flame. The second combustion mode (CASE B 
mode) corresponded to the regimes where flameholding became 
useless and burning can self -develop . The parametric calculations 
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were done to find location of the boundary which separates 
operational conditions corresponded to different combustion modes. 
Results demonstrated that transition boundary located at very 
high level of the incoming air static temperature (Ta=2200-2500K) . 
Its decreasing to the air temperature level which could be of any 
practical interest (Ta=1700-1900K) required significant increasing 
of the H 2 temperature (Tf ~ 900 - 1000K) . Even for scramjet 
operational conditions, estimations demonstrated that transition 
from CASE A to CASE B combustion mode could be expected at very 
high flight Mach numbers (M^-14-16) . Thus, this finding 

demonstrates that, now, performed calculations have mostly 
academic interest for the study of the flamelet solutions features 
in near-critical regimes rather than any practical applicability. 
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APPENDIX A 


ACCOUNT OF THE FLAMELET EQUATIONS (1.2) 


The outlined account of the flamelet model equations was 
proposed by V. Kuznetsov in [21] . 

Let us consider the instantaneous mass conservation equation 
for the reactive specie 

p -|^+p (UV) Ca=VpDVCa+pRa (A.l) 

where p is density , U is flow velocity in the laboratory frame; 
Ca is the reactive specie mass fraction; D is its diffusivity; Ra 
is chemical production term. 

Let us introduce the mixture fraction conservation equation as 
dz 

P +P (UV) z=VpDVz (A. 2) 


Consider the stoichiometric surface z=z =l/(l+St) where St is 

S 

the stoichiometric coefficient. The velocity of this surface V in 
the laboratory frame is given by the relation: 

v= Oz/at-vz) / (VzVz ) 

Using eq.(A.2) one obtains: 

V= (U-Vz) Vz/ (VzVz) - i (VpDVz) Vz/ (VzVz) 

Let us define quantity: 

W= ( (V-U) Vz) / (VzVz) 1/2 

which is normal to the surface z = z g component of the vector V in 
the frame moving with a media. It is given by equation: 

W = ^ (VpDVz) / (VzVz) 1/2 


The quantity W has the Kolmogorov scaling since it depends only on 
gradients of the scalar field. Therefore one has estimate: 

.. , . 1/4 

W ~ (ev) 

where c is the mean dissipation, v is the kinematic viscosity. 

Let us choose some point x=x on the surface z = z and some 

0 s 

time instant t=t Q and adopt new frame moving with the velocity 
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W= ( (V-U) Vz) / (VzVz ) 1/2 

The eqs . (A.l) could be written as 

p +p (WV) Ca=v P DVC « + P Ra 

Let us make now some estimates assuming the reaction zone 
thickness 6 to be small . One has 


VpDVCa ~ pD|^ ~ 

^ an K an „2 

O 

where n is normal direction to the surface z=z 

C 

1/4 


One has also 


pWVCa 


pWCa 

5 


pCa (ev) 
<5 


Since D~i> one can conclude that the molecular diffusion term 
VpDVCa is much larger than the convective term pWVCa if 
5<< T)= {V /£) 

where r) is the Kolmogorov scale. 

The same kind of reasoning could be applied to the 
non- stationary term p . Hence eqs. (A.l) can be reduced to a 
form: 


a=l , ,J 


(A. 3) 


VpDVCa + pRa = 0 

Small thickness of the reaction zone allows to apply two 
additional simplifications: i) neglect the variation of pD inside 
this zone; ii) consider zone as locally planar and to neglect the 
derivatives along the stoichiometric surface. Therefore eqs. (A. 3) 
are simplified still further: 


D Rot = 0 

dn 2 


oc = X , . . . , J 


where n is coordinate which is normal to the surface z=z . 

s 

Finally the mixture fraction is introduced as an independent 
variable instead of n using the relation z=z g + (g^- ] n. As the 

^ ' s 

result one has the flamelet model equations in a form (1.2 
Sec .1.1: 


of 


N 


S 


d 2 Ca 

dz 2 


+ Ra 


0 


a=l ( . . . J 
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APPENDIX B 


DESCRIPTION OF THE SOLVERS USED IN THE RESEARCH WORK 


B . I . FLAMELET EQUATIONS SOLVER FLSLV 


The numerical solution of the flamelet model equations (1.5) 
was based on time relaxation method. The finite difference 
approximation was based on the central difference approximation 
for the second derivatives and the linearization of the chemical 
production term over reaction species mass fractions. As the 
result the finite difference system of equations in delta form 
was as follows (afterwards subscripts denote grid points along z 
and superscripts denote levels along t) : 


SC" 


N a 


2SC 


25C J 


2SC J 


i+l 


(Az +Az ) Az 

i + l i i 


Az Az 

i+l l 


( Az +Az ) Az 

i i+l i+l 


= N 


2C 


j -i 

i -1 


(Az +Az )Az 

i + l i i 


j-1 


2C 


Az Az 

i + l i 


2C 


J - i 
i + l 


(Az +Az ) Az 

i i+l i+l 


w J V 

i 


aw 

-» 

ac 


j-i ^ 
SC J 


where x J =(t j -t J *) is pseudotime step; Az =(z -z ) is 

-> -4 -4 i i i-i 

C ) is vector of the 


in z-direction; SC J = (C J 


(B.l) 


grid step 
dependent 


variable increment; OW/3C) is Jacoby matrix. 

The obtained system of the finite difference equations was 
block- tridiagonal . At each pseudo-time step system (B.l) was 
solved using matrix sweep method where the elementary block 
matrixes were inverted using Gauss reduction. When the reactive 
species mass fractions were obtained at j-th pseudotime level, 
the energy conservation equation was used to calculate 

temperature distribution. The nonlinear algebraic equation for 
temperature was solved using simple iterations. The iterations were 
controlled by the relative error ER= max mod [ (hcaic - hz ) /hz ] 
where maxima was got through all grid points z i (i = l,...,l). 
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Static enthalpy h^ was calculated using obtained at the 

ca c i 

current time level values of reactive species mass fractions and 
the static temperature obtained in the current iteration on 
temperature. Static enthalpy value h^ was calculated directly as 
function of z i from the dependence between enthalpy and mixture 
fraction of (1.5) . The iterations were stopped when ERclO -3 . 

The increments SC control was applied for choice of the optimal 

pseudo- time step r and to fasten the convergence to steady state 

solution. The algorithm was as follows. The norm for SC was 

calculated at each time step as c J = max max{mod(SC ) J /(C ) J_1 ) 

oc i 1 a i a i 1 

where maxima was got through all the species (a=l,...,ii) and all 
grid points z (i=l,...,l). Two bounds e=0.5 and e=l.O were 

i i u 

fixed. If e was greater than e u then time step was diminished by 
a factor of 2.2 and calculations were repeated for the same time 
level again. If e J <e u then time step remained the same and 

calculations of the next time level was performed with the same 
time step. If c ^ was smaller than then time step was increased 
by a factor of 2 and calculations for the next time level were 
performed. Maximum value of time step was restricted by z =10. 

max 

The time relaxation was stopped when the value of e J became lower 
then £ q =10 6 ■ The obtained solution was considered as converged. 


B . 2 . 2D PARABOLIZED NAV I ER- STOKES EQUATIONS SOLVER SUPNEF 

The governing PNS equations were solved by marching method 
using numerical code SUPNEF, which was previously developed in 
CIAM for the steady supersonic chemically reacting flows [51] . 
The code is based on the explicit finite difference method which 
is the generalization of the well known steady analogy of 
Godunov scheme [38] for the steady supersonic flows. The method 
is based on the conservation laws for the elementary 
computational cell. The modification [52] , [53] of Godunov scheme 
providing for Euler equations the higher order accuracy both 
in lateral and longitudinal directions is realized. In current 
case the predictor-corrector version [54] is used. The finite 
difference scheme is of the second order accuracy on regular 
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uniform grids and retains approximation on arbitrary 
nonuniform grids. The higher order accuracy scheme is monotone 
due to using principle of minimal increments of functions in the 
cell [55] . This principle is analogy to well-known minmod 
limitator [56] . The solution of two supersonic flows 

interaction problem generalized to the multispecies flow case is 
used at the corrector step to approximate the convective fluxes 
on the cell boundaries. The viscous tenses are approximated 
explicitly by the values in the previous cross-section using 
central differences modified to include nonuniform grids. 
Formally such viscous terms approximation decreases the order of 
approximation up to the first one. However this scheme 
allows to exclude the main approximation errors associated with 
the convective inviscid terms. The explicit scheme has the 
limitations on step of integration in longitudinal direction. 
However it is not a serious limitation for the pure supersonic 
jet-like flow. 

The implicit approximation was used for the chemical source 
terms in species mass conservation equations. The value of 

W a (C,p,T) in conservation equation for a-specie was represented in 
a following form: 

W a (C,p,T)=-f* (C,p,T)C a +f‘ ( C , p , T) 

so, that f* and f^ were positive functions. This expansion 
provides the solution monotonicity of the model equation 

dt = -£ c + f • 

The implicit approximation of the chemical source terms 
required to use the iterative method for the solution of the 
finite difference equations system for each cell. The iterative 
procedure included two loops: the global iterations which were 

necessary to find the gasdynamics parameters at known 
concentrations (at the finish of internal iterations) and the 
internal iterations for concentrations calculations to take into 
account the influence of temperature and density variations 
obtained in the previous global iteration. In the global iteration 
step the gasdynamics subsystem is solved at the known mass 
fractions. Internal iterations are carried out to define the mass 
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fractions from the species concentrations subsystem for the known 
(from previous global iteration) corrected temperature and 
density. The internal iterations for mass fractions were 
performed using Gauss-Seidel iterations. Internal iterations on 
concentrations are performed up to the maximum absolute accuracy 
less than 10 5 . Then the new correction of gasdynamics parameters 
is provided (the following global iteration) and so on up to the 
convergence of iterations both for thermodynamics parameters and 
mass fractions. The iterations convergence for gasdynamics 
parameters is controlled by the relative increment of temperature. 
This value is chosen to be less than 10” 5 . If the number of 
iterations becomes too large the longitudinal step of integration 
is twice decreased. 

B . 3 . 2D STEADY NAVIER- STOKES EQUATIONS SOLVER FNAS2D 

The system of averaged NS equations was solved numerically 
using modified version of FNAS2D code developed at CIAM [42] . It 
is based on the time relaxation procedure for the modified version 
of Godunov's scheme providing second order accuracy of the steady 
state solution on regular near-uniform grid and approximation on 
arbitrary nonuniform grid. The method is realized on the grid 
constructed in physical plane without transformation to the 
computational plane. 

The system of conservation equations is written for each cell 
with implicit approximation of convective and viscous fluxes 
through boundaries using time increments of main dependent 

variables vector where C was C T =(p,U, V, P, v , K, z,a 2 ) for FL 

approach calculations and C T = (P, U, V, P, , {C a , a=l , . . . , J} ) for QL 
approach calculations. 

The arbitrary discontinuity breakdown problem (which is used 
in original Godunov scheme for convective fluxes approximation on 
the cell boundaries using parameters in cell centers) is used in 
modified form to include them into implicit scheme. The Riemann 
problem solution for the new time level is considered in linear 
approach. It is supposed that the configuration of main 
discontinuities realized after breakdown on the new time level is 
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identical to this one on the old time level. Therefore the 
identical relations exist between the parameters on the cell 
boundaries with the parameters in the cells centers before 
breakdown on the old time level from one side and between 
parameters increments in time on the boundaries and in the cells 
centers from the other. If the pressures ratio on the main shock 
and/or on the expansion wave is higher than 1.2 the iterative 
procedure is used to obtain the accurate solution of the nonlinear 
Riemann problem on the old time level. These corrected values are 
used for the convective fluxes approximation on the cell 
boundaries on the old time level without correction for the new 
time level . 

To provide the higher order accuracy for the steady state 
solutions the piecewise linear parameters distributions within 
cell are assumed and the parameters in the middle point of each 
cell boundary edge are estimated for adjacent cells on the old 
time level to provide the initial data for the Riemann problem. 
These parameters are evaluated with the aid of minimal derivatives 
principle [55] modified to arbitrary nonregular grids in [57] . The 
minimal derivatives or minimal space increments principle provides 
the monotonicity condition [58] of the higher order accuracy 
scheme. As to the parameters increments in time, their values from 
the cell centers are used as the initial data for the Riemann 
problem (it is acceptable for time relaxation procedure without 
steady state solution accuracy loosing) . 

The viscous tenses and diffusion fluxes on the cell boundaries 
on the old time level are approximated with the aid of central 
differences generalization for arbitrary grids. The contribution 
of viscous terms into finite difference operator on the new time 
level is taken into account approximately because the contribution 
of some grid points is dropped out to provide the five diagonal 
structure of the algebraic system of equations for time increments 

of the dependent variables §C- 

The implicit approximation is used for the chemical sources in 
reactive species conservation equations (for QL calculations 
only) . 

The resulting system of algebraic equations for the main 
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(B . 2 ) 


parameters increments had the following structure: 

-> -> -> -> -> 

A 6C +B SC +C SC +D SC +E SC =R 

i,j i, J s i+i,j i,j s i,j i,J S,J-1 i,j S.j+l i,J 


here i and j are indexes for the current grid cell center; A, B, 
C, D, E are matrices with size (8,8) in the case of FL approach 
and (5+J , 5+J) in the case of QL approach, J is total number of 

reactive species; R is the residual of the steady state equations. 
The system (B.2) was solved in the following manner: 

i. Jn the case of FL approach calculations - First the equations for 

~ 2 

v , z, cr and K (turbulence subsystem) were solved independently 
one from another on new time level. Further the effective heat 
capacities ratio F and effective heat of formation Q on the new 
time level were approximated using flamelet library and were 
averaged. Then the gasdynamics subsystem (I I. 7) was solved taking 
into account the increment of the enthalpy in time estimated with 
the aid of known increments in time of f and Q. 

ii. In the case of QL approach calculations - The full coupling in 
the implicit part of the difference operator between gasdynamics, 
concentrations and turbulence model subsystems was realized. This 
techniques provided much more better convergence characteristics 
of the implicit scheme than the sequential solution of the 
aforementioned subsystems reported in our interim Report [27] . 

Algebraic systems of finite difference equiations were solved 
using point Gauss-Seidel method. Sequential downstream and 
upstream sweeping in Gauss-Seidel method were performed for each 
time level. The number of internal iterations ( in Gauss-Seidel 
method) on each global iteration in time was chosen depending on 
the global residual level and Courant number. 
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Appendix C. 


Table C.l. 


Reaction rates constants for hydrogen oxidation chemistry based on data by 

Miller and Bowman [32], 


Reaction rate constant is presented in the form K = AT 8 exp(- — — ) . 


Units used are Kelvins, seconds, mols and centimeters. 


RT 



Reaction 

Forward rate 

lg A 

B 

-E/R 

1 

H2 + 02 = OH + OH 

13.23 

0.00 

-24044.0 

2 

OH + H2 = H20 + H 

9.07 

1.30 

-1824.7 

3 

0 + OH = 02 + H 

14.60 

-0.50 

0.0 

4 

0 + H2 = OH + H 

4.70 

2.67 

-3165.3 

5 

H + 02 + M = H02 + M 

17.56 

-0.72 

0.0 

6 

OH + H02 = H20 + 02 

12.88 

0.00 

0.0 

7 

H + H02 = OH + OH 

14.15 

0.00 

-540.0 

8 

O + H02 = 02 + OH 

14.15 

0.00 

-540.0 

9 

OH + OH = 0 + H20 

8.78 

1.30 

0.0 

10 

H + H + M = H2 + M 

18.00 

-1.00 

0.0 

11 

H + OH + M = H20 + M 

22.20 

-2.00 

0.0 

12 

H + 0 + M = OH + M 

16.79 

-0.60 

0.0 

13 

0 + 0 + M = 02 + M 

13.28 

0.00 

899.8 

14 

H + H02 = H2 + 02 

13.10 

0.00 

0.0 

15 

H02 + H02 = H202 + 02 

12.30 

0.00 

0.0 

16 

H202 + M = OH + OH + M 

17.11 

0.00 

-22896.7 

17 

H202 + H = H02 + H2 

12.20 

0.00 

-1912.2 

18 

H202 + OH = H20 + H02 

13.00 

0.00 

-905.8 

19 

N + NO = N2 + 0 

12.51 

0.30 

0.0 

20 

N + 02 = NO + O 

9.81 

1.00 

-3160.2 

21 

N + OH = NO + H 

13.58 

0.00 

0.0 
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Table C.2 


Reaction rates constants for hydrogen oxidation chemistry based on data by Warnatz [34], 

£ 

Reaction rate constant is presented in the form K = AT B exp (-—~ ). 

RT 

Units used are Kelvins, seconds, mols and centimeters. 



Reaction 

Forward rate 

lg A 

B 

-E/R 

1 

H02 + H = OH + OH 

14.18 


-505.1 

2 

H02 + H = H2 + 02 

13.40 


-348.8 

3 

H02 + O = OH + 02 

13.30 

0,00 

0.0 

4 

H02 + OH = H20 + 02 



0.0 

5 

H02 + H02 = H202 + 02 

12.30 

0.00 

0.0 

6 

OH + OH + M = H202 + M 

22.11 


0.0 

7 

H + H202 = H2 + H02 

12.23 

0.00 

-1888.3 

8 

H + H202 = H20 + OH 


0.00 

-1804.0 

9 

0 + H202 = OH + H02 

13.45 


-3223.4 

10 

OH + H202 = H20 + H02 

12.85 

0.00 

-721.6 

11 

H + 02 = 0 + OH 

17.08 


-8311.0 

12 

H + H + M = H2 + M 

17.81 

-1.00 

0.0 

13 

H + OH + M = H20 + M 

22.93 

-2.00 

0.0 

14 

H + 02 + M = H02 + M 

17.85 

-0.80 


15 

OH + H2 = H20 + H 

8.00 

1.60 

-1659.8 

16 

OH + OH = O + H20 

9.18 

1.14 

0.0 

17 

02 + H2 = OH + OH 

13.74 


-29106.5 

18 

02 + M = 0 + 0 + M 

18.26 


-59415.5 

19 

O + H2 = H + OH 

7.18 



20 

N2 + 0 = NO + N 

14.04 


-37496.7 

21 

N + 02 = NO + 0 

9.81 


-3160.2 

22 

N + OH = NO + H 

13.58 


0.0 
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Table C.3 
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Fig. 1 . Schematic of the flamelet approach. 
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Fig. 3. Scheme and conditions of Beach et.al. coaxial experiment. 
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Fig.4. Flamelet library generation. 

Residual norm e and time step behaviour. N— 15 sec'! 
Solution for N s = 10 sec' 1 is used as initial conditions 
for time relaxation. 
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Fig. 6. Estimation of calculations accuracy for H and OH mass fractions 
based on 'three grids' extrapolation. Calculations for the case N— 100 sec' 1 . 

I is the number of grid points. 

mean root square interpolation. 
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Fig. 10. Flamelet model calculations based on: 
— Miller — Bowman detailed kinetics; 

Wamatz detailed kinetics. 

N = 100 sec' 1 ; Beach test case. 


0.5 

0.4 

0.3 

0.2 

j 

b 0.1 
0.0 

1.0 


H20 mass fraction 



H20 mass fraction 



0.00 0.20 0.40 0.60 0.80 1.00 

mixture fraction z 


Fig. 1 la. Sensitivity of flamelet solution to the exponent p value in correlation (1.8) 

Water mass fraction. 

N®=50 sec' 1 . Beach test case. 
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Fig. lib. Sensitivity of flamelet solution to the exponent p value in correlation (1.8) 

OH mass fraction. 

INr= 50 sec" 1 . Beach test case. 






Fig. 12. Computational domain together with the adaptive grids used in calculations 

(100 cells in y- direction) 
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Fig. 14a. Pitot pressure P t profile in the initial cross section (x/dj = 0.33). 
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Beach test case. 
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Fig. 14b. Centerline distribution of the hydrogen mean mass fraction. 

FL approach calculations 

^ experimental data 


Beach test case. 
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Fig. 15a. Influence of the cells number on the results of FL computations 

Beach test case, x/dj =8.26 
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Fig. 15b. Influence of the cells number on the results of QL computations. 

Beach test case, x/dj =8.26 
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Fig. 16. H 2 O mass fraction contours 
Beach Test Case 
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Fig. 17. Calculated distribution of N® vs distance from the nozzle. 

Beach test case. 
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Fig. 18. Mach number contours 
Beach Test Case 
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Fig. 19. C >2 mass fraction behavior in the vicinity of the self- ignition point 
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Fig. 20. Turbulent mixing characteristics contours 
FL approach 
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Fig.21a. Comparison of FL and QL predictions together with experimental data. 

FL approach; 

QL approach; 

experimental data of Beach et al. 
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Fig. 21b. Comparison of FL and QL predictions together with experimental data. 

FL approach; 

QL approach; 

experimental data of Beach et al. 
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Fig.21c. Comparison of FL and QL predictions together with experimental data. 

FL approach; 

QL approach; 

^ experimental data of Beach et al. 
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Fig. 2 Id. Comparison of FL and QL predictions together with experimental data. 

FL approach; 

QL approach; 

^ experimental data of Beach et al. 

x/d, = 27.9 
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Fig. 22. Scheme and conditions of Burrows and Kurkov experiment. 
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Fig. 24a. Flamelet library. Role of pressure variation. 
Water mass fraction. 

P s = 0.08 MPa 
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P s = 0.12 MPa 
Burrows-Kurkov test case. 
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Fig. 24b. Flamelet library. Role of pressure variation. 
"Effective" heat of formation. 

P s = 0.08 MPa 
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P s = 0.12 MPa 
Burrows-Kurkov test case. 
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Fig. 24c. Flamelet library. Role of pressure variation. 
"Effective" heat capacities ratio. 

P s = 0.08 MPa 
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P s = 0.12 MPa 
Burrows-Kurkov test case. 
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Fig. 25. Sensitivity of flamelet solution to detailed chemistry model. 
Flamelet model calculations based on: 

Miller - Bowman detailed kinetics 

Wamatz detailed kinetics 

N s = 30 sec'! P s =0.1 MPa. 

Burrows-Kurkov test case. 
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Fig. 26b. Sensitivity of flamelet solution to the exponent p variation. 

Water mass fraction. 

N s = 30 sec' 1 P s =0.1 MPa. 

Burrows-Kurkov test case. 
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Fig.26d. Sensitivity of flamelet solution to the exponent p variation 
"Effective" heat capacities ratio T 
N s = 30 sec" 1 . P s =0.1 MPa. 

Burrows-Kurkov test case. 



Fig.27. Computational domain together with grid (90x100 cells) 
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Fig. 28. Profiles at the inlet boundary. 
Burrows — Kurkov test case. 
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Fig. 30. Velocity profile resolution in the vicinity of the lower duct wall. 

Exit cross section (x=0.356m). 

FL approach. 

O - cell centers location. 

Burrows- Kurkov test case. 
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Fig.31. Convergence history. 

FL approach 
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Fig. 32a. Influence of discretization on result of 
FL predictions for species mole fractions. 
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Fig. 32b. Influence of discretization on result of 
QL predictions for species mole fractions. 
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Fig. 33. Water- Vapor mass fraction contours 
obtained in calculations 


Burrows - Kurkov test case 
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Fig. 34. Comparison of FL and QL predictions together with experimental data. 

x=0.356 m. 

FL approach 

QL approach 
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Burrows-Kurkov test case. 



Fig.35. Comparison of FL and QL predictions together with experimental data 

x=0.356 m. 
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Fig. 36. Relative role of the chemical nonequilibriumness/mixture fraction fluctuations 

in FL approach predictions budget. 

Flamelet averaged solution. 

Flamelet instantaneous solution. 

Equilibrium chemistry limit 
Experimental data 

Burrows- Kurkov test case. 
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Fig. 37. Results sensitivity to variation of pdf shape 

Airy/Gaussian P(z) model 
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Fig. 39. Results sensitivity to the turbulent 
fluctuations intensity level INT=cr/ 
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Fig.40. Computational expences. 



Fig.41. Comparative estimation of computational expences. 

Beach test case. 


a) Beach experiment b) Burrows-Kurkov experiment 
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Fig. 43. Accuracy of approximate relation for effective oxidation rate. 

0<N S <100 sec ; P = 0.1 MPa 
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Fig. 44. Accuracy of a C* calculations based on approximation (III. 17). 
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Fig.45a. Accuracy of simplified PEq method. 

accurate FL solution 

simplified PEq method 


Beach test conditions. 
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Fig. 45b. Accuracy of simplified PEq method. 

accurate FL solution 

simplified PEq method 


Beach test conditions. 
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Fig.45c. Accuracy of simplified PEq method. 

accurate FL solution 

simplified PEq method 
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Fig.45d. Accuracy of simplified PEq method. 
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Fig. 46a. Accuracy of simplified PEq method. 

accurate FL solution 

simplified PEq method 

Burrows — Kurkov test conditions. 
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Fig.46b. Accuracy of simplified PEq method. 

accurate FL solution 

simplified PEq method 
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Fig. 46c. Accuracy of simplified PEq method. 

accurate FL solution 
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Burrows — Kurkov test conditions. 
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Fig. 47. Different combustion modes predicted by flamelet model. 
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Fig.48a. Critical values of scalar dissipation vs air 
and fuel temperatures. 

P = 0. 1 MPa 
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Fig. 48b. Pressure influence on and N^! 
Tf=900 K 
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Fig.49. Boundary between combustion modes vs regime parameters. 
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Fig.50. Comparison with typical scramjet operational conditions. 
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